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Abstract 

This thesis describes a study of a cavity QED microlaser in which many atoms are 
present simultaneously and atom-cavity interaction is well-defined. The microlaser is 
found to display multiple thresholds analogous to first-order phase transitions of the 
cavity field. Hysteresis is observed as a function of atom-cavity detuning and number 
of atoms. Data is compared with a rate equation model and fully quantized treatment 
based on micromaser theory. Good agreement between theory and experiment is 
found when the cavity is resonant with atoms of the most probable velocity, but long 
lifetimes of metastable states preclude the observation of true steady-state transition 
points. For nonzero atom-cavity detuning the microlaser displays broadenings and 
shifts which are not yet well-understood. 

Quantum trajectory simulations are performed to investigate many-atom and fi- 
nite transit time effects in the microlaser. We show that over a wide range of parame- 
ters the many-atom microlaser scales with the single-atom theory, with a perturbation 
in the photon statistics due to cavity decay during the atom transit time. 
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Chapter 1 



Introduction 

1.1 Cavity QED and the microlaser /micromaser 

One of the simplest models for light-matter interaction consists of a single two-level 
atom interacting with a single electromagnetic field mode of a cavity. If the atom- 
cavity coupling is stronger than the atomic or photon decay rates, the irreversible 
decay of an excited atom is replaced by an oscillatory exchange of energy between 
the atom and cavity, in a manner characteristic of two coupled oscillators. 

In the 1950's this model was the subject of considerable theoretical study, primar- 
ily in connection with the then-recently invented maser pQ . Experimental realization 
of such a strongly coupled system would not occur until several decades later, however, 
with advances in Rydberg states and superconducting microwave cavities. 

In the early 1980's the micromaser, or single-atom maser, was invented by Walther 
and collaborators [2j. In this experiment, individual Rydberg atoms pass through and 
emit photons into a superconducting microwave cavity with extremely high quality 
factor. The micromaser formed a close approximation to the single-atom, single- 
mode ideal. Subsequent experiments showed that the micromaser exhibited a variety 
of interesting and unusual phenomena, including trap states, near-number states, 
bistability, and nonclassical statistics [H]-[7]. 

The recent development of multilayer dielectric mirrors with very low loss has 
made it possible to conduct similar experiments at optical wavelengths. In 1994 
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Kyungwon An and Michael Feld developed the microlaser, the optical analogue of 
the micromaser [HI E] • This device is also known as the cavity QED microlaser in 
order to distinguish it from other microlasers, for example those using semiconductor 
microresonators. 

This thesis describes the first experiments exploring the microlaser's truly unusual 
properties. In particular, we show that in contrast to conventional lasers which display 
a single laser threshold, the microlaser exhibits multiple thresholds. 

1.2 Microlaser concept 

The basic scheme of the microlaser is illustrated in Fig. 11-11 A beam of two-level 
atoms pass through a resonant high-finesse optical cavity. Before entering the cavity 
mode each atom is excited to its upper state via a laser 7r-pulse. The atom interacts 
with the mode for a time t mt , during which it may emit a photon into the cavity. 

The key characteristic of the microlaser is not only the strong atom-cavity coupling 
but that different atoms experience almost the same interaction with the cavity: both 
the atom-cavity coupling strength g and the atom-cavity interaction time t int are well- 
defined. It will be shown in chapter El that the resulting emission probability can be 
written as 



where n is the number of photons present in the cavity prior to the entrance of the 
atom. The unusual properties of the microlaser are primarily due to this sinusoidal 
dependence of the atomic emission probability on the cavity photon number. 

1.3 Previous work 
1.3.1 Micromaser 

Reviews of the micromaser can be found in a number of references in the bibliography. 
Most relevant for this experiment is the study by Benson et al ^T] of dynamic behavior 



P< 



emit 
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Figure 1-1: Schematic of cavity QED microlaser. A: atomic beam, M: mirrors, C: 
cavity mode, P: pump field, 8: cavity tilt angle (cf. [TU]). 

near the micromaser field phase transitions. Spontaneous jumps and slow transitions 
between metastable states are observed; hysteresis is seen as interaction times are 
varied. The results were explained by quantum Monte Carlo simulations and a Fokker- 
Planck description of the micromaser dynamics. The presence of small stray fields 
prevented a full quantitative agreement between theory and experiment from being 
reached; the calculated transition rates were generally much slower than what was 
observed experimentally. 

1.3.2 Microlaser vs. micromaser 

We comment here on how the microlaser differs from its microwave counterpart. 
Many atoms required 

The threshold atom number as derived in Chapter El is 

A^th = (1.2) 
For our current microlaser with supersonic beam we have T c « 180 kHz, t- m t ~ 
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0.1/xsec, and g ~ 193kHz, giving iV th ps 10. 

This is to be compared with the micromaser for which typical figures are T c « 
30 Hz, ti n t ~ 47/xsec, and g m 40 kHz (values from ^1]), so threshold atom number 
N t h ~ 3 x 10~ 4 , almost five orders of magnitude lower. 

One consequence of the many-atom nature is that the trap state resonances, cor- 
responding to situations in which atoms perform an integer number of entire Rabi 
oscillations, cannot observed; nor can the associated near-number states of the micro- 
maser jjj. However, other interesting effects appear from the many-atom case. The 
many atom nature of the current microlaser is treated in detail in Chapter |3] 

No thermal photons 

At finite temperature every field mode has a nonzero photon occupancy n tn given by 
the Planck formula. For the micromaser operating at ~ 20GHz, blackbody photons 
nth is as high as 0.15 even at 0.5 K and have a significant effect on the micromaser, as 
discussed in many micromaser articles. In the microlaser at 300K, n tn ~ 10~ 25 and 
can be ignored. 

Constant transit time 

Micromaser experiments are typically performed with a fixed atomic injection rate and 
variable interaction times. This roughly corresponds to varying the pump parameter 
9 as described in micromaser theory (see Chapter^. 

Use of a supersonic beam of the microlaser prevents any significant variation in the 
transit times. Instead, the parameter most easily varied is the number of intracavity 
atoms N. As N is varied, both the expectation photon number n ex and the pump 
parameter 9 change. 

Direct detection of photons 

Optical frequency photons can be easily outcoupled out of a cavity, and have suffi- 
ciently high energy to be directly individually via a photomultiplier tube or avalanche 
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photodiode. This is in contrast to the micromaser experiments, in which maser oscil- 
lation must be inferred via observation of the atom states. Direct detection of photons 
opens up new possibilities for photon statistics or laser linewidth experiments. 

Atom-cavity coupling 

The centimeter-scale wavelengths used in the micromaser allow the atomic beam to 
largely avoid coupling nonuniformities due to mode structure. In the micromaser, 
standing-wave coupling variation requires us to employ a traveling-wave interaction 
scheme (see Chapter 0]). In addition, the small dimensions of the cavity mode and re- 
quirement of uniform coupling place strict requirements on the position and geometry 
of the atomic beam coupled to it. 

1.3.3 Optical cavity QED 

At optical wavelengths the experiments most similar to this one have been spec- 
troscopy and probe transmission experiments of high-finesse optical cavities contain- 
ing a small number of atoms. 

Heinzen and Feld [121 E] performed the first observation of enhanced and su- 
pressed spontaneous emission and level shifts in an optical cavity. 

The vacuum normal- mode splitting of atoms in a cavity was observed in [T^j . 

Oscillatory transmission and nonclassical photon statistics of a probed cavity was 
reported by the Kimble group (T^j. Optical bistability was also observed. 

To our knowledge, all these optical experiments employed atomic excited states 
with lifetimes ~ns much shorter than atom-cavity transmit times. They were per- 
formed in passive (non-laser) cavities in which a probe beam was incident. In addi- 
tion, the atom-cavity couplings were not well-determined due to standing-wave mode 
profiles (cf. chapter HJ). 
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1.3.4 Lasers and optical bistability 

Bistability in the microlaser bears some similarity to "classical" optical bistability 
(see, for example, [T6.]) in active and passive cavities. Multimode laser bistabilities 
due to mode competition in gas lasers have been studied since the earliest days of 
the laser jTTj. Current-feedback bistability has been observed in some semiconductor 
laser systems, e.g. [TS] . 

1.3.5 Phase transition analogies 

A detailed description of a laser is necessarily statistical in nature, due to the presence 
of randomness in many forms and large numbers of degrees of freedom. It is therefore 
natural to draw analogies between laser behavior and other problems in statistical 
physics. In particular, for some time physicists have drawn analogies between a laser 
near threshold and matter near a second-order phase transition. This comparison was 
first described in detail by DeGiorgio and Scully in 1970 PHJESI], although mention 
of this analogy date back considerably earlier (cf. references in ibid.). 

In the present work we investigate jumps in the microlaser field which can be 
shown to be analogous to first-order (discontinuous) phase transitions [21. To our 
knowledge no similar behavior has been observed in any other laser. 

1.3.6 Many-atom dynamics in cavity QED 

The field of collective effects in radiative properties goes back at least to Dicke |22j . 
who showed that a collection of iV two-level radiators, localized within a distance 
smaller than their radiative wavelength, may emit photons at a collective rate that 
scales with N 2 instead of the usual N. 

Many-atom effects in the micromaser have been studied by a number of authors 
(ES]-j2S])- It is well-known that the micromaser's trap state resonances are easily 
destroyed by two- atom events ([2E] and [2*7j). 

A study by Kolobov and Haake j2H| of two-atom events predicts a change in 
average photon number as well as location of the first-order thresholds, relative to the 
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single-atom theory. On the other hand, Elk [21], considers a micromaser injected with 
multiple atoms, equally spaced in time, and finds a scaling behavior very similar to 
the single-atom theory. D'Ariano [22] considered a micromaser pumped with clusters 
of up to several hundred atoms at once, an interesting model but one that cannot be 
considered realistic in the present context. 

The issue of many-atom effects is central to this thesis and is considered in de- 
tail in Chapter El We find a result similar to that of [21], but we have a differing 
interpretation of an observed extra variance. 
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Chapter 2 

Theory of the microlaser 



We begin by considering the ideal model of a single atom in a cavity. The microlaser 
is introduced by an semiclassical rate equation model which describes the basic fea- 
tures of the microlaser in an intuitive way. We then give a fully quantized treatment 
following the density matrix approach [28J. A Fokker-Planck analysis from [21] is re- 
viewed which reproduces the rate equation solutions and also describes the transitions 
between different solutions in terms of an effective potential. 

2.1 Atom in a cavity 

The Hamiltonian describing the interaction of a radiation field with a single-electron 
atom is 



where 1~La and Jip are the atom and radiation field energies, r is the position of the 
electron, and E is the electric field operator: 



U = H A + H F - er ■ E 





k 



H A = -huJ a CTz 
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(2.3) 



E = 5>k£k(cik + 4) (2.4) 



k 

Here the sum is over all allowed wave vectors k and 



is a function of mode volume 14 which for a near-planar cavity is given by 

V = ^w 2 m L (2.6) 

where L is the cavity mirror separation and w m is the Gaussian mode waist at the 
center of the cavity. 

The atom-field interaction is given by 

n AF = %kV+ak + 4<r_) (2.7) 

k 

where 

gj> = e JMiL^ (2 . 8) 

For a two-level atom, we have in the rotating-wave approximation (i.e. neglecting 
non-energy-conserving processes) 

H = Y^ Mc4°k + \^a<y z + %k(o-+a k + 4°"-) ( 2 - 9 ) 
k 1 k 

where a^, are photon creation and destruction operators, a + , a + are atom 
creation and destruction operators, huj a is the atom level spacing. 

For a single mode cavity of frequency u c we have 



H = Tiw c (va H — Uoj a a z + hg(a + a + aV_) (2-10) 
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To go to the interaction picture we write 



H =H + H 1 

H = hu c a^a + \Tiu) a a z (2-H) 
H 1 = hg(cr + a + aV_) 

The Hamiltonian in the interaction picture is 

V = e iHot/h Hie -iH t/h = hg ( a+ae i*t + a i a _ e -i^ (212) 

where A = oj c — u a is the atom-cavity detuning. 

ih^\1>) = V\1>). (2.13) 

The atom-cavity wavefunction can be written 

00 

\m) = EW*)|e,n> +c g>n (t)\g,n)} (2.14) 

n=0 

where |e, n) corresponds to a state with an excited atom with n photons in the cav- 
ity, etc. Clearly the Hamiltonian mixes only the states |e, n) and \g,n + \). Therefore 
the Schrodinger equation can be written as the coupled equations 



c e ,n = -igVn + le zAt c gjn+1 , (2.15) 

c g ,n+i = -igVn + le~ tAt c Cin (2.16) 
A general solution is given by 

C e ,n(t) = {ce,n(0) C0s(fi n t/2) - *— sin(fi n i/2) (2.17) 



2 ' !!X/I! V -c 9iri+1 (0)sin(^/2))e^/ 2 
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c g,n+l(t) — \Cg,n+l(Q) 



iA 

cos(Q n t/2) + — sin(fi n t/2) 



(2.18) 



2igy/n + 1 



a. 



c e>n (0) sin(n n t/2) e^ 2 



where 



n n = V A2 + % 2 (n + l) 



(2.19) 



If the atom is initially in the excited state, then c e>n (0) = c n (0) and c 9i „ + i(0) = 0, 



and 



c e , n {t) = c„(0) 



iA 

cos(tt n t/2) - — sin(a*/2) 



3 iAt/2 



(2.20) 



Cg,n+l(^) = -Cn(0) 



sm(aV 2 ) 

i in. 



-iAt/2 



(2.21) 



Consider the case of perfect resonance A = and initial condition of excited atom 
with exactly Uq cavity photons. The solution is then given by 



5 ,no(*) = cos(\/no + 1 0*) 



Cg, no +i(t) = sin(Vno + l gt) 



(2.22) 
(2.23) 



with all other c e>ri (t), c ff) „(i) equal to zero. 

2.2 Rate equation for microlaser 

Eq uat ions 12 . 221 and 12 . 231 describe resonant quantized Rabi oscillation. Note that oscil- 
lation proceeds even for an initially empty cavity (no = 0) (vacuum Rabi oscillation). 

In the microlaser, atoms interact with the cavity for an interaction time tint. The 
emission probability for an atom is given by the squared ground state amplitude after 
its interaction with the cavity: 



^emit = |c 9 ,n +l(^int)| = SU1 (y/no + 1 fft) 



(2.24) 
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where for simplicity we suppose that the cavity initially contained exactly no photons. 

Let us now write a rate equation for the cavity photon number n: the time rate 
of change of n is equal to to the difference between gain and loss terms: 



^ = G -L = ^- sin 2 (VnTT gt-m) -T c n = (2.25) 

at tint 

which is set equal to zero at steady-state. Here N is the number of atoms in the 
cavity, which can be defined by N/t- m t = r, the injection rate. Loss is equal to the 
photon number multiplied by the cavity loss rate T c . 

The rate equation can be solved graphically as in Fig. 12-11 In general there exists 
more than one solution. 

A solution n' is called stable if the system responds to a small perturbation in n 
in a manner that tends to restore the original solution. 

d 

Stability condition: — (G - L)L> < (2.26) 
on 

In terms of Fig. 12-11 a solution is stable if and only if the slope of the loss line 
exceeds that of the gain line at their intersection. 

The threshold atom number N ta is given by the point at which the gain and loss 
terms have equal slope, assuming small argument of the sine function: 

= tL (2 - 27) 

Solutions to !2.25l as the injection rate r is varied are plotted in Fig |2-2I . with gt int = 
0.1. The number of solutions increases indefinitely with increased r. Solutions are 
asymptotic to constant-n lines corresponding to integer number of Rabi oscillations 



y/n + 1 gt- m t = twk, m integer. 

The rate equation analysis shows how the atoms' sinusoidal emission probability 
leads to bistability and multistability. However, its derivation was less than rigorous 
and assumes that the cavity photon number is well-determined. In addition, it gives 
no information about the relative stability of different solutions, or the cavity photon 
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1000 2000 3000 4000 5000 6000 7000 8000 
photon number n 

Figure 2-1: Dimensionless gain (solid line) and loss (dashed line) in the semiclassical 
rate equation analysis for experimental g, i mt , and N = 1000. Closed circles represent 
stable solutions, open circles unstable solutions. 
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Figure 2-2: Multiple solutions of microlaser rate equation. Dotted lines represent 
photon numbers corresponding to an integer number of complete Rabi oscillations 
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statistics. 

For a more complete theory we turn to the quantum theory developed for the 
micromaser by Filipowicz, Meystre, Scully and others 



2.3 Quantum theory of microlaser 

We follow the density matrix approach [2E|- This theory in based on the following 
assumptions: 

1. Excited atoms are injected at a rate low enough so that at most one atom at a 
time is present in the cavity. 

2. Excited and "ground" state lifetimes are much longer than t int . 

3. Cavity damping during interaction time is small: r c t int <C 1 

2.3.1 Density matrix equation of motion 

Gain term 

Consider the change in the cavity field density matrix due to the injection at time to 
of an excited atom, which then interacts with the cavity for time t int : 



The state p nn '{to + tint) m &y be found by determining the atom-field evolution from 
time to to to + tint and then tracing over atomic variables: 



Note that density operators with four indices represent combined atom-cavity sys- 
tems; those with two represent the field only. A coarse-grained derivative for the gain 




(2.28) 




(2.29) 
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contribution to the field density matrix can now be written 
/ dpnn' \ 



dt 



— r [P en;en' (to + tint) + Pgn;gn'(to + ^int) ~~ Pen; en' (to)} (2.30) 
gain 



We now replace to with the arbitrary continuous time t. The time-evolved density 
matrices p en;en '(t + t int ) and p gn - gn >(t + t int ) are determined from the solutions for the 
excited state amplitudes as in section 1 of this chapter. We have 



C e ,n(t + tmt) = C en (t) COs(V U + 1 gt- mt ) 

c g ,n+i(t + t- mt ) = Cg tn+1 (t) sin(Vn + 1 gt int ) 
which gives for the density matrices 

Pen\en'(t + ^int) = c e,n(t + tj n t)c* n '(t + tint) = 



Pnn'(t) COs( V / n + 1 #t m t) COs( V / ^' + 1 fitfint) 
Pgn;gn'(t + ^int ) = c g,n(t + ^int) c ^ )?1 ' (t + ^int ) = 

Pn-l,n-l'(t) COs(v/n gt int ) COs(v / ^ 7 gt int ) 

We then have the reduced field density matrix equation of motion: 

(^Sr) ■ = _r t 1 ~ COS(y™ + 1 9tint) C0S( V / ^' + l9tint)]pnn' + 
\ aL / gam 

r sm(y/ngt int ) sin( v / n 7 fi'tint)pr i -i,n'-i- 



f2.31^ 



(2.32) 



(2.33) 



Note that we have made no assumptions about the photon number distribution 



Loss term 

The contribution due to cavity loss is 



1 ' "' 1 (a) ap — 2apa) + pa) a) (2.34) 



dt L 2 

/ loss 



for zero thermal photons. This is a standard result from the theory of oscillators 
coupled to a reservoir (c.f. (2H| chapter 8.) 
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Combining 12.331 and I2.34[ the total density matrix equation of motion is 



Pnri — a nn'Pnn' + &n-l,n'-lPn-l,n'-l + c n+l,n'+lPn+l,n'+l (2.35) 



where 



'[1 - COS<y™ + 1 S^int) COS(Vn' + lfltfint)] ~ if (™ + 



6 nn / = r sin(- v /n^ int ) sin(\/n 7 5 , t 



int J 



(2.36) 



If we restrict ourselves to diagonal matrix elements this master equation reduces to 



Pnn QnnPnn "n—l,n—lPn—l,n—l Cn+l,n+lPn+l,n+l 



(2.37) 



where 



a nn = ~r sin 2 (v / rTTT gt- mt ) - T c n 
b nn = r sin 2 {^/ngt int ) 



(2.38) 



2.3.2 Steady-state photon statistics 

Setting p n (t) = p nn (t) = gives the equation 

{r sin 2 (Vw + 1 fltfint) + T c n} p(n) + r sin^y^iatMn - 1) + 
r c (n + l)p(n + 1) = 

which leads to the following recurrence relations: 

{r sin 2 (Vn + 1 gt int ) + T c ri}p(n) = T c (n + l)p(n + 1) (2.40) 

r sm 2 {^/ngt mt )p{n — 1) = T c np(n) (2-41) 
33 



We obtain the steady-state photon distribution of the microlaser: 

f\ (n\ TT r sin 2 (Vkgt int ) fo An . 

p(n) = p{0) [[ — (2.42) 

k=l 1 c K 

where p(0) is determined from normalization. 

5>(n) = 1 (2.43) 

n 

It is convenient to define two quantities, the expectation photon number 

N cx = r/r c (2.44) 

equal to the number of atom injections per cavity decay time, and the normalized 
interaction time or pump parameter 



9 = \fN ex gt int (2.45) 

which is the Rabi phase if N ex photons are present. Linearization of the rate equation 
for small n and setting gain equal to loss gives 9 = 1 as the laser threshold condition. 

Mean photon number 

In Fig. 12-31 the photon number has been plotted as function of pump parameter 9 for 
N ex = 10. Note the dips due to integer number of Rabi oscillations (trap states) and 
a "jump" near 9 = 2tt which becomes sharper as iV ex becomes large. 

Photon statistics 

The width of the photon number distribution can be characterized by the Mandel Q 
parameter, defined as 

Q - ^^7^ - 1 (2-46) 
(n) 

The Q parameter is equal to zero for Poisson statistics (for which the variance (n 2 ) — 
(n) 2 = (n), greater than zero for super- Poisson statistics (e.g. thermal light) and less 
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than zero for sub-Poisson statistics. 

In Fig. 12-41 the Mandel Q parameter is plotted as function of pump parameter 9 
for N cx = 10. The cavity photon statistics are sub-Poisson over wide ranges of 9. 

To compare with the rate equation result we now calculate the results from 
Eq. 12.421 as the number of cavity atoms is varied, with fixed gt- 1T &. 

Fig. 12-51 shows the average photon number from Eq. 12.421 together with the rate 
equation result. We see that the fully quantized theory agrees with one of the rate 
equation stable solutions except during "jumps" corresponding to transitions between 
the solutions. The quantum theory essentially "chooses" one of the rate equation 
solutions. 

During the transitions, which become increasingly narrow as N e x is increased, the 
photon number distribution has two peaks corresponding to the two stable solutions. 

It can be shown that the rate equation solutions correspond to within ±1 of a local 
maxima of the p(n) distributions from the quantum theory. Suppose a local maximum 
in p{n) exists at n — n ; this means that p(n ) > p(n Q ± 1) and the quantity inside 
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Figure 2-4: Micromaser theory result for Q, with N ex = 10. Q < indicates sub- 
Poisson statistics. 
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Figure 2-5: Rate equation solutions (dotted lines) compared with quantum microlaser 
theory (solid line). 
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Figure 2-6: Q parameter according to microlaser quantum theory. Statistics are 
strongly sub-Poissonian for almost every point over threshold. Spikes are due to 
bimodal photon distributions during phase transitions. 



the product of Eq. 12.421 passes through 1 between k = n and k = n + 1 : 



rsin (Jn^gtint) 1 rsin (y/n + lgt mtJ 



T c n T c (n + 1) 

Setting 



> 1 > ; v / u ' i " mv (2.47) 



rsin (y/np + lgt int ) = . , 

r c (n + i) 1 ' ) 



gives the original rate equation Eq. 12.251 



2.4 Modification for realistic parameters 

We describe modifications to the rate equation and quantum theory in order to ac- 
comodate the atomic velocity distribution, nonuniform atom-cavity coupling, atom- 
cavity detuning, and imperfect pumping. 
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2.4.1 Averaging of gain function 



The framework for all these modifications is the same: we calculate the effect of the 
relevant broadening on the "gain function" f3(n) which for an ideal (monovelocity 
etc.) microlaser is given by 

#>(n) = sinVn+1 gt int ). (2.49) 
For example, the effect of a velocity distribution f v (v) is the average 

_ POO 

p(n) = / sm 2 (V^+Tgt int (v))f v (v)dv (2.50) 
Jo 

The steady-state photon statistics are then given by 

p(»)=p(°)n^jr ( 2 - 51 ) 

k=l L c K 

Note that we essentially replace a process in which different atoms experience dif- 
ferent interactions, with a process in which all atoms experience the same averaged 
interaction. This may seem implausible until we realize that the end result 12.511 rep- 
resents an average itself. Our averaging procedure follows from the assumption that 
the random variables associated with the velocity, etc. are independent (cf. |2*T]). 

The following expression is used to give the average of the gain function over 
velocity distribution, coupling variation, and detuning: 



(3(n) = / B(n,v,A(v)) f v (v) f g (g) dv dg (2.52) 

J 

where A(t>) = kv6 and B(n, v, A) is the gain function including detuning. It might be 
expected that B(n,v, A) is given by the well-known solution involving the generalized 
Rabi frequency 

B(n, v, A) = Af ^ n ^ 1 } 9 A9 sin 2 U(n + l)g 2 + A 2 t int l (wrong) (2.53) 



4(n + l)g 2 + A 2 
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However, this expression is only valid for constant fields. In general B(n,v,A) must 
be calculated numerically and its form will depend on the shape of the coupling 
profile g(x = vt), not just its area. The function B(n, v, A) is found by integrating 
the Bloch equations ( Eqns. l2~T5l and l2*.16j) for a initially excited atom in a field with 
detuning A and Rabi frequency y/n + lg(t) during a simulated transit through the 
cavity. Integration of B(n,v, A) for every calculation would be very time-consuming; 
therefore we tabulate this function on a grid and construct an interpolating function 
to be called on by our theory calculations. 

2.4.2 Imperfect inversion 

The adiabatic pumping process described in Chapter 0] is not perfect: for the data 
shown in the Chapter the excited state probability was about 80%, measured by a 
fluorescence experiment. This includes a small contribution from atomic decay during 
the pump process and the transit time between the pump field and cavity mode. 

Let p ee and p gg be the upper and lower state probabilities, respectively. We claim 
that for N ex very large the effect of nonzero p gg is simply to reduce the effective 
number of atoms to 

N cS = N( P ee-p gg ). (2.54) 

assuming p ee > p gg . That is, instead of the number of atoms we may speak of the 
inversion. This idea is a familiar one in laser physics. 

To support this claim we note that atoms initially in the ground state have a 
probability of removing (absorbing) a photon from the field of 

-Pabsorb = sin 2 (v/n gt mt ) (2.55) 

Intuitively it is clear that the case p ee = p gg = 1/2 would produce essentially zero 
photons, since the initially lower state atoms absorb photons in the same manner as 
the upper state atoms emit them. 

39 



G,L 

1 r 




1000 2000 3000 4000 5000 6000 7000 8000 



n 

Figure 2-7: Dimensionless gain (solid line) and loss (dashed line) in the semiclassical 
rate equation analysis, including effects of velocity distribution, nonuniform coupling, 
and detuning, with N = 1000. Closed circles represent stable solutions, open circles 
unstable solutions. 

In the rate equation picture, the net emission rate per atom is then modified to 

-Pemit = Pee SlD. 2 (y/n + 1 fltffait) - Pgg ^(y/n gt int ) « (p ee - p gg ) Sin 2 (\/n+ 1 gt int ) 

(2.56) 

with corrections on the order of 1/2-^/n. Corrections may affect the initial threshold 
region somewhat but for most of the range of interest in our experiments, n is very 
large. 



Realistic rate equation and quantum theory 

Fig. 12-81 gives the equivalent of Fig. 12-51 after realistic effects, appropriate to the data 
presented in Chapter have been included. 

The average photon number can also be plotted as a function of detuning. (Fig. 12- if 
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Figure 2-8: Photon number vs. atom number, rate equation (dashed line) and micro- 
laser theory (solid line) including effects of velocity distribution, nonuniform coupling, 
and detuning. 
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Figure 2-9: Q parameter according to microlaser quantum theory, including real- 
istic effects as in Fig. 12-81 Spikes are due to bimodal photon distributions during 
transitions. 
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Figure 2-10: Rate-equation solutions for N e g = 1000. Second and third branches 
appear as closed curves. Upper halves of closed curves represent stable solutions. 
Small gaps are artifacts from algorithm used to find solutions. 

2.5 Fokker-Planck analysis 

In [21] a Fokker-Planck equation for the cavity photon number distribution p(n, t) is 
derived: 



0_ 

Of 



d 1 d 2 

p(n,t) = -—Q(n)p(n,t) + - — [G(n)p(n,t)\ 



(2.57) 



where for zero thermal photon number 



Q(n) = R sin 2 {y/ngt mt ) - T c n 



(2.5* 



G(n) = R sin 2 (y/ngt- mt ) + T, 



(2.59) 



and R is the injection rate. To generalize to our averaged case we replace the 
sin 2 ( v /n^t int ) terms with j3(n) from 12.521 (The replacement of n by n + 1 is in- 
significant on the photon number scale we will consider.) Equation 12.571 was shown 
to reproduce the average photon number and photon statistics of micromaser theory 
as long as (i) the pump parameter 6 {nN cyi ) l l 2 and (ii) n ^> 1. Both of these 
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Figure 2-11: Solutions to rate equation with realistic parameters, as a function of 
detuning, N cS = 600. Dotted line, quantum theory shows jumps in detuning curve. 



conditions are easily met in our case. 

We now rescale time to the cavity lifetime and photon number to expectation 
photon number: r = tT c , v = n/N ex . Then Eq. 12.571 becomes 

d d 1 d^ 1 

di P ^ t] = ~Yv q{y)v{y ' r) + 2^ (Z7)p(Z7 ' r)] (2 ' 60) 

with 

q(y) = f3(vNex) - v (2.61) 

g {y) = (3(uNex) + v (2.62) 

noting that A^ ex = R/T c . 

The stationary solution to the Fokker-Planck equation is then 

p(u) = ~~~7~t exp ( 2N CX r dv (2-63) 

where C is a normalization constant. For large iV ex the photon number distribution 
will accumulate in the global maximum of the exponent of Eq. I2.63| or equivalently 
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the global minimum of an effective potential defined by 



V(u) = -f 




q{u) 



(2.64) 



The local minima of V(u) are the zeros of q{v), which are the solutions of 



= v 



(2.65) 



which is equivalent to the rate equation derived in 12.251 

We now have an intuitive picture for the transitions between stable points: they 
occur as the global minimum of the effective potential V(v) changes from one point 
to another. This is a variation on the Landau theory of phase transitions, with y/v 



sharp in a similar manner to the sharpness of phase transitions for a large number of 
particles. 

We note that the initial threshold is of second-order: as the injection rate is in- 
creased there is a continuous increase of the global minimum in the effective potential 
from v = 0. 

2.5.1 Metastability and hysteresis 

In general the microlaser does not automatically find the global minimum but must 
fluctuate to it from out of a local minimum. Therefore hysteresis is expected if 
parameters are varied faster than the transition rate between metastable and stable 
states. 

To give an example, let us consider the rate equation solutions as detuning is 
varied. As shown in Fig. 12-111 the first branch appears as an inverted parabola- 
type shape centered around zero detuning (relative to resonance of most probable 
velocity atoms); higher branches first appear as closed curves, and with increased 
atom number, join with the first branch solution. (See Fig. 12-12)1 . Suppose that 
metastable lifetimes are very long and spontaneous transitions do not occur. Then 



as the order parameter. Note that as N ( 



cx 



oo the transitions become infinitely 
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Figure 2-12: Rate equation solutions (realistic parameters) as a function of detuning, 
iV = 2500. Arrows show transitions that could occur due to disappearance of solutions 
as detuning is changed. 



as detuning is modulated we would expect a consistent hysteresis pattern following 
the arrows in Fig. !2-l2|) . These transitions occur when parameters change so that a 
local minimum in the potential disappears entirely, and the system is forced to the 
next potential minimum it encounters (not necessarily the global minimum). 

Similar features would be encountered for photon number as a function of atom 
number. 



2.5.2 Transition rates 



The tunneling rate from a metastable solution to a global minimum is estimated via 
the Kramers analysis [21] to be, in units of cavity linewidth, 



W 



1 

27 



\q{VM)— ( — r 



1/2 



exp{-2iV ex [V>M) - V{v m )]} 



(2.66) 
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Figure 2-13: Tunneling from local to global minimum in the effective potential V(v), 
plotted for realistic microlaser parameters with N c g = 600. 

where v m is the metastable local minima point and Vm is the local maximum over 
which the system must tunnel to reach a global minimum (or possibly another, lower 
energy, metastable local minimum). (Fig. 12-13*)) 
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Chapter 3 

Numerical simulations of 
many-atom microlaser 

We present results from quantum trajectory simulations of a microlaser/micromaser 
device in which more than one atom at a time may be present and atom transit times 
are not required to be much smaller than the cavity decay time. For parameters 
in which the system is well removed from the influence of trap states the average 
photon number is in good agreement with an appropriately scaled single-atom, weak- 
decay micromaser theory. The photon number variance, however, shows an increase 
proportional to N/N ex = r cav t int for parameters in which the weak-decay theory 
predicts sub-Poisson statistics.. This result is similar to that found by Elk in the 
case of regular atomic injection j21]. By considering the case N ex = 1 we show that 
this broadening of the photon number distribution is not a many-atom or correlation 
effect, but is due to cavity decay during the atom's transit through the cavity. 

3.1 Introduction 

The microlaser and micromaser differ from conventional lasers and masers not only 
in that they use single atoms as the gain medium, but in the controlled nature of the 
atom-cavity interaction. Uniformity of the coupling strength between atom and cavity 
[TO] and interaction time £i nt of the atom through the cavity, in addition to the long 
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lifetime of the atom upper level state compared to the transit time, create an unusual 
situation in which different atoms undergo nearly the same interaction with the cavity. 
That is, interaction is almost independent of decay and other random effects. This is 
in contrast to descriptions of a conventional laser [20] in which ensemble averages are 
performed to account for the lifetimes of the excited atoms, inhomogeneous effects, 
etc. 

It is natural to ask which properties of the microlaser/micromaser persist even 
when the restriction of a single atom is lifted. Dynamics based on single-atom quan- 
tum Rabi oscillations are not necessarily expected to hold; in addition, correlation 
effects may cause unknown differences from single-atom theory. 

A related problem is to determine the effect of abandoning an assumption made 
in analytical treatments of the microlaser/micromaser j2I] that the cavity decay time 
is much longer than the interaction time, i.e. r cav t int << 1, where T cav is the cavity 
linewidth. In this paper we will call this the condition of weak decay. 

The present study was motivated by the realization that with the current ex- 
perimental parameters, more than one atom at a time must be present for laser 
oscillation to occur. In going from the standing wave to the tilted atomic beam, 
uniform-coupling, traveling-wave interaction ^U] the coupling strength g is decreased 
by a factor of 2 relative to the peak standing-wave value. (This factor of 2 was incor- 
rectly reported as a factor of y/2 in JU|.) The threshold number of atoms is given by 
EH: 

Nfh = Fcav ^ 4 

y tint 

This is to be compared with iV t h ~ 1 for the original (standing wave) microlaser 
experiment [8 j if g is taken at an antinode. Therefore the study of multiple atom 
effects is necessary for understanding the microlaser with traveling-wave interaction. 

Effects relating to finite r cav t int are also of interest in the microlaser, since with 
present parameters [TQJ we have r cav t int ~ 0.1. 
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3.2 Previous approaches 



The effect of 2-atom events was studied in [73] in the context of micromaser trap 
states. It was shown that collective events destroy trap state resonances with great 
efficiency. 

More general treatments of multi-atom behavior have involved strong simplifica- 
tions to the atom injection scheme. 

D'Ariano studied a micromaser system pumped by clusters of up to N = 100 
atoms which all enter and leave the cavity at the same time j2S]- It was shown that 
for iV < N ex the system behaves similarly to the one-atom maser, and for N > N ex 
the system exhibits multiple thresholds. 

Elk [21] considered an injection scheme in which atom arrivals are equally spaced 
in time. Results for small 9 and N < N ex showed a nearly constant average photon 
number but a linear increase in photon number distribution width o with N/N ex . 
However, it is difficult to determine whether the results are related to regular injection 
of atoms, which in general causes a reduction of noise (i.e. smaller a). We note that 
the results in the case iV = 1 do not agree with single-atom micromaser theory |21j . 

In this paper we study the case of realistic, Poisson atomic injection, and interpret 
the results based on finite transit-time rather than many-atom effects. 

3.3 Quantum trajectory analysis 

To simulate the microlaser with random (Poisson) injection of atoms, we used a 
quantum trajectory algorithm j221 5 [S3]- This technique was previously used to analyze 
the threshold-like transition, and the possibility of many-atom effects, in the original 
standing- wave microlaser in [3U|. This paper contains for a detailed description of the 
extended quantum trajectory technique also employed in the present study. 

Due to limitations of computer memory and processing time, for each simulation 
run it was necessary to set a maximum value iV max for the number of atoms that 
could be present in the cavity at a given time. In each case iV max was adjusted such 
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that the "overflow" probability was 1% or less, ensuring that the statistics are very 
close to Poissonian. 

In single-atom micromaser theory [21], the behavior of the system is determined 
by three parameters: (i) the expected photon number iV ex = r/r cav where r is the 
atom injection rate; (ii) the pump parameter 9 = y/N ex gt int where g is the atom-field 
coupling constant; and (iii) the thermal photon number n b . 

It is then shown that, with the assumptions of single-atom events and r cav t int << 1 
(weak decay), [21] the steady-state photon number distribution in the cavity is given 
by 



n b ™ n 

Pn =P0 



l + n b 



n 



k=l 



l N ex s in 2 (VA; + lgt mt ) 



(3-1) 



rifjk 

where po is determined by normalization. In this paper Eq. 13.11 will be referred to as 
the result from micromaser theory, although it is also applicable to the microlaser. 

In this study we set thermal photon number n b = 

We performed runs of 10,000 - 100,000 atom injections, for varying values of N ex , 
9 and number of intracavity atoms iV = rt int . We scaled g and r c as necessary to 
obtain the correct values for iV ex and 9 as N was varied. At the beginning of each 
run, at least 2,000 atom injections were neglected in order to allow the system to 
reach steady state before data collection began. 

Most simulations were performed with N ex = 10. This value was chosen for several 
reasons. First, a relatively low photon number allows larger numbers of atoms to be 
considered. Second, effects that scale as N/N ex are more apparent for smaller N ex . 
Finally, we observed that simulations reached their steady-state distributions much 
more rapidly for lower iV ex . 

For simplicity, and to allow comparison with analytic theory, we regard g and ti nt 
to be perfectly well-defined (i.e. uniform coupling and monovelocity atoms). Our 
QTS program has also allowed us to simulate the effects of various types of variation 
and broadening in a realistic system. Results of these investigations will be reported 
elsewhere. 
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theta 

Figure 3-1: Two-atom events destroy the trap state resonances 

The simulation was written in C (see Appendix) and all calculations were per- 
formed on a 400 MHz Pentium II workstation running Linux. The calculations pre- 
sented here required roughly 2000 hours of CPU time in total. 

3.4 Results 

3.4.1 Single atom limit; Trap states 

Fig. 13-11 shows the average photon number as a function of 9, for N ex = 10. Note 
the deep depressions due to trap states resonances in which atoms perform integral 
numbers of Rabi oscillations. 

With N = 0.1 the QTS results for average photon number (n) show significantly 
smaller modulation due to trap states than the result of Eq. 13.11 although the proba- 
bility of the cavity containing 2 or more atoms at any given time is ~ 2%. It is only 
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Figure 3-2: Quantum trajectory results for average photon number 

with iV = 0.01 that QTS results give good agreement with micromaser theory. This 
illustrates the sensitivity of the trap state resonances to 2-atom events as described 
in [731. 

The photon number distributions in the low- A/" limit were verified to be in excellent 
agreement with Eq. 13.11 further confirming that the simulation was working properly. 



3.4.2 Broadening of photon number distribution 

Figs. 13-21 and 13-31 show values of average photon number n and Mandel parameter 
Q = ((n 2 ) — (n) 2 )/ (n) — 1 as a function of 9 for a range of atom numbers N, with 
N ex = 10. Note that for larger values of 9, Q is large and may be difficult to interpret 
due to the presence of 2 or more peaks in the photon number distribution. 

For 9 < 5, before the trap state resonances occur, the results for average photon 
number remain remarkably close to the prediction of micromaser theory. The values 
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Figure 3-3: Quantum trajectory results for Mandel Q parameter 
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Figure 3-4: Universal curve showing an increase in Q which scales with N/N ex 

for Q, however, show a significant increase with N in the range 2 < 8 < 5. This 
interval corresponds to the region in which according to micromaser theory the system 
displays sub-Poisson statistics (Q < 0). 

Similar increases in Q are observed in a small interval near 8 = 8, between two 
trap states. Simulations with N ex =30 and 50, in which trap state resonances are 
more narrow, showed similar behavior over wider intervals between trap states. 

The parameter Q was calculated as a function of N/N ex , with 8 = 3, for iV ex = 5, 
10, 15, 30. The increase of Q is seen to be linear with N/N ex ; for N ex > 5 we find 
(Fig l3-4|) that all points may be approximately described by 

Q = Q (N ex ,8) + a(8)(N/N cx ) (3.2) 

where Qq is the value of the Mandel parameter from micromaser theory and a(8 = 
3) w 0.80 ± .05. 
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Figure 3-5: Slope a (9) of the increase in Q with N/N ex , as a function of pump 
parameter 
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Figure 3-6: Increase of variance with N/N ex for iV ex = 1 



Fig. 1123 shows a(8) calculated via a least-squares fit with iV cx = 10 and < iV < 2. 
Error bars indicate the degree of linearity. It is seen that a exhibits oscillation with 
8 and is positive for 9 where the single-atom, weak-decay theory predicts sub-Poisson 
statistics (Q < 0). 



3.4.3 Many-atom vs. finite transit time effects 

An increase of variance proportional to N/N ex was previously found by Elk, who 
interpreted it as being due to many-atom interactions [24 . 

We note, however, that N/N ex is equal to simply ti nt T c , a parameter describing 
cavity decay during the transit time, and is not directly related to the number of 
atoms. The extra variance observed in this study and by Elk can therefore be seen 
as a consequence not of multiple atoms, but of cavity decay during an atom's transit 
time. For example, we may have N ^> 1 but virtually no extra variance if N ex ^> N. 
More importantly, there are examples in which iV <C 1 but excess variance still occurs. 
Fig. 13-61 shows Q as a function of N for N ex = 1,6 = 2. A linear increase in variance 
is seen even for iV = 0.1. 
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Figure 3-7: Dashed line: upper bound of Q parameter, including finite decay time 
effect. Solid line: Q parameter from microlaser theory without without finite decay 
time effects, as in Fig. 12-91 

3.5 Predictions for experiment 

To apply our quantum trajectory results to our microlaser system, we make a seem- 
ingly precarious assumption that simulations for N < 5 atoms represent the behavior 
for up to iV ~ 1000. If true then the predictions for the experiment are as follows. 

No significant deviation from the single-atom theory was found for the few-atoms 
microlaser; therefore we would expect the cavity photon number to follow the pre- 
diction of the microlaser with realistic parameters described in Sec. 12.41 as plotted in 
Fig. El 

For the photon statistics, we have in the microlaser N/N ex w 0.1, and therefore 
the perturbation to photon statistics is fairly small; it would not differ strongly from 
that given in Fig. 12-91 

An upper-bound prediction for the Q parameter is given in Fig. 13-71 
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3.6 Discussion 



Our calculations have clarified the distinction between many-atom effects and ef- 
fects due to cavity decay during transit-time. Many-atom events destroy the trap 
state resonances and therefore have strong effects for parameters for which these res- 
onances are prominent. Away from these resonances, the system shows remarkable 
agreement with single-atom micromaser theory. However, for values of 9 for which 
micromaser theory predicts sub-Poisson statistics, cavity decay during transit time 
cause an increase in variance relative to micromaser theory, by an amount propor- 
tional to N/N cx = t int T c . 

Lack of collective effects away from trap states may be a consequence of the 
noncorrelated Rabi phases of atoms which have entered the cavity at different times. 
Increase of variance represents an absolute decrease in the negative differential gain in 
the atom-cavity system. Future work will include development of models to explain 
the observed effects quantitatively. 
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Chapter 4 



Apparatus 



and methods of 



microlaser 



experiment 



In this chapter we describe the components of the microlaser and techniques used in 
its study. In particular we describe how a well-defined interaction between atoms and 
cavity can be achieved. 

4.1 Two- level atom 
4.1.1 Barium energy levels 

Figure EFT] shows the levels of 138 Ba relevant to the microlaser experiment. The mi- 
crolaser uses the 1 So *-> 3 Pi transition of wavelength A = 791.1 nm and linewidth ap- 
proximately 50 kHz. The excited state decays to the ground state and two metastable 
D states with a branching ratio of 0.43:0.41:0.16. We use this intercombination line 
due to its relatively long lifetime T" 1 = 1.3 /isec . 

The most frequently studied transition in 138 Bais 1 So «-> x Pi at 553 nm. It has 
a linewidth of about 19 MHz. The fluorescence spectrum of a naturally occuring 
barium sample is shown, with isotope shifts, in Fig. 14-21 

The 553 nm transition was used to measure atomic density in the cavity, determine 
pump beam efficiency, and measure the atomic velocity distribution. 



59 




Figure 4-1: Relevant energy levels for 138 Ba. 



Natural occuring barium, of which 71.7% is 138 Ba, is relatively inexpensive and 
readily available. We purchased the metal in rod form, packed in argon, from Alfa 
Aesar. 

4.2 Optical resonator 

The microlaser cavity mirrors consist of superpolished fused silica substrates coated 
with a multilayer, ion-sputtered dielectric coating. The substrates and coatings were 
purchased from Research Electro-Optics (Boulder, Colorado). The coating consists 
of 45 alternating layers of Si0 2 and Ta 2 5 . 

4.2.1 Characteristics 

The microlaser cavity consists of two mirrors in a symmetrical, near-planar Fabry- 
Perot configuration. The mirror separation, measured via transverse mode spacing, 
is 1.1 mm. The radius of curvature of the coated surfaces is 10 cm. 

The mirrors have an overall diameter of 7?mm and length of ? mm. The actual 
reflective surfaces are 45-degree beveled to a diameter of 3 mm in order to facilitate 
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Figure 4-2: Isotope spectrum of the 1 Sq *-> x Pi transition at 553 nm. 
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alignment of the tilted atomic beam and cleaning of the coated surfaces. 

The mirrors are glued into stainless steel holders (Fig. 14- 3|) designed to hold the 
mirrors without birefringence-inducing stress. The cavity is aligned outside the cham- 
ber, then the holders are glued on the inside of a cylindrical piezo transducer (PZT) 
used to modulate the cavity mirror spacing. 
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Figure 4-4: Microlaser resonator assembled with PZT 



The microlaser experiment is performed with the fundamental (TEMoo) mode of 
the cavity (it is also possible to use higher-order transverse modes). The mode waist 
is given by 

w = (^pj « 41.1/xm (4.1) 

where r is the mirror radius of curvature, L is the cavity spacing, and A is the 
wavelength. 

The cavity mode volume is given by 

V = \Lw\ (4.2) 
and the peak atom-cavity coupling for the standing-wave case is 



9o= M2^ = 385 kHz (43) 



The atom-cavity interaction time for a particular atom is defined as 

= ^ (4.4) 

V 
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Figure 4-5: Support system for cavity assembly in chamber. 
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Figure 4-6: Cavity transmission of an earlier cavity (from 0) showing TEMOO and 
higher order transverse modes. Transverse mode spacing in the current cavity is 7.11 
GHz. 

where v is the atom's velocity. The a/7t factor accounts for the Gaussian shape of the 
cavity mode: the integral of g with time for a moving atom is the same as an atom of 
the same velocity traversing a "top hat" mode of constant coupling g for time t- mt : 

/+oo 
g exp 
-00 

The cavity spacing L was determined via measurement of the transverse mode 
spacing. In a near-planar cavity, 

Az/ t = L > r (4.6) 

2tttJLR/2 

where higher-order corrections can be shown to be ~ 10~ 3 for our cavity [S2 ch. 19]. 

We determined Au t by measuring the difference in PZT voltages between the fun- 
damental and first transverse mode resonances. The scan was calibrated by measuring 
the PZT voltage shift between two sidebands of the 24.67 MHz FM modulation used 
for laser frequency locking (see Sec. X). We find the PZT scan calibration is 171.6 
MHz/V and verified that it does not vary measurably over the PZT voltage range in 
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Figure 4-7: Schematic of optics for microlaser experiment. To the right of the cavity 
a mirror on a "flipper" mount is used to switch between cavity fluorescence detection 
for velocity distributions and detection of cavity transmission. 

question. The transverse mode spacing was measured to be 



Au t = 7.11 GHz 



(4.7) 



giving a cavity length 



L = 0.90 mm. 



(4.8) 



4.2.2 Cavity PZT 

The PZT controlling the cavity separation is a tube (Vernitron PZT-5A) of height 
0.75" and inner diameter 0.75", and wall thickness 0.125". Four holes are drilled in 
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log signal 




Figure 4-8: Log of ringdown signal (a.u.) vs. time. Line shows least-squares linear 
fit. Here T cav = 1.02 //sec and finesse F = 1.0 x 10 6 . 

the cylinder wall to allow passage of the atom and pump beams. 

Isolation mechanisms are used to isolate the cavity from vibrations. Viton o-rings 
hold the PZT inside a heavy brass housing (Fig. 14- 3|) . The housing is placed in 
a supporting structure by fittings made of teflon (Fig. I4-5J) : the support in turn is 
secured inside the vacuum chamber with teflon screws. 



4.2.3 Finesse measurement by ringdown technique 

The finesse of a Fabry- Perot cavity is defined as the ratio of free spectral range (FSR) 
to the FWHM of a cavity transmission peak, or equivalently. 

F = ^ n ^L- (4.9) 
l-R l-R y 1 

where R is the geometric average of the intensity reflectivity of the two mirrors, 
i.e. R = y/RiR~2 and we assume R very close to unity. 

To measure finesse we scan the cavity PZT through resonance with a probe laser 
beam. The cavity decay time is measured by fitting a curve to the exponential 
decay of the cavity field (Fig. I4-8|) . The signal shows a sinusoidal modulation due 
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Figure 4-9: Figure 4. Anomalous cavity transmission. Scan speeds: (i) 0.6, (ii) 4.8 
GHz/s, compared with theory; (iii) unmodified cavity lineshape. (b) Frequency shift 
vs. inverse cavity scan speed. The best fit is for A = 0.2 ppm. 



to interference between the probe beam and the Doppler-shifted cavity field. This is 
described further in Ref. 



4.2.4 Absorption measurement via thermally-induced opti- 
cal bistability 

The finesse is given by tt/(1 — R), with R « 1 the mirror reflectivity and 1 — R the 
total loss, equal to the sum of three loss parameters, the absorption A, scatter S 
and transmission T. R, hence the total loss, can be measured by cavity ringdown 
[HE]- However, determining how the loss is distributed among A, S and T is difficult. 
In particular, few techniques exist for measuring very small absorption coefficients. 
Methods such as photothermal deflection |37j are capable of measuring A down to 1 
ppm. 

Recently we have developed a new technique to measure A at the sub-ppm level 
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[3H1 EH] • This technique relies on the intense circulating field created when a resonant 
probe beam is incident on a high finesse cavity. Any absorption of the intense internal 
field causes heating in the coating layers, and the subsequent thermal expansion of 
the mirror substrate reduces the mirror separation distance and leads to a thermally- 
induced optical bistability of the empty resonator. This optical bistability causes 
distortion in the measured cavity lineshape, as observed by slowly scanning the cavity 
length. The lineshape exhibits an extended resonance with a slow rise and rapid 
drop (Fig. 4), the duration of the resonance being directly related to the absorption 
coefficient A [38j. 

We have characterized all three parameters and established that A « 0.2 ppm and 
S ~ 2.5 ppm, thus establishing that scattering is the predominant loss mechanism in 
our supercavity mirrors 
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Figure 4-10: Microlaser output as a function of cavity-atom detuning, for varying 
atomic beam tilt angles 9. From 

4.2.5 Uniform coupling by traveling-wave interaction 

In the original microlaser experiment, the atom beam is incident normal to the cavity 
axis and atoms interact with the standing wave of the Fabry-Perot cavity. Since the 
coupling g is proportional to the local field amplitude, which varies as sin 2 (fcz), atoms 
experience a wide range of coupling strengths. 

We eliminated this problem by introducing a small tilt 9 <C 1 between the atom 
beam direction and the normal to the cavity axis. Due to the Doppler shift, each 
atom then experiences not one standing wave but two traveling waves of frequencies 
io ± kv9. If 9 is sufficiently large, the atom will be resonant with only one traveling 
wave component at a time. This eliminates the fast spatial variation due to the 
standing wave. It should be noted that in the traveling-wave configuration the peak 
coupling constant, g , is reduced by a factor of 2 compared to the peak standing- wave 
value (see Appendix |DJ). (This factor of 2 was incorrectly reported as \/2 in [TUj.) 
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4.3 Atomic beam 



A typical atom beam apparatus consists of a container in which the sample is heated 
in vacuum to form a vapor; a small nozzle allows vapor atoms to escape. Apertures 
are placed downstream to collimate the beam. 

We have two main requirements for our barium atomic beam: high beam density 
and narrow velocity distribution in the ground state. 

4.3.1 Velocity selection: mechanical, optical, or supersonic 

An atomic beam with narrow velocity spread is needed to provide a well-defined 
atom-cavity interaction time (i.e. uniform transit time), and to ensure that all atoms 
traversing the pump laser field prior to entering the cavity are excited to a state of 
complete inversion (tt pulse condition). 

The "classical" technique of mechanical velocity selection (via e.g. Fizeau wheel) 
has several serious drawbacks: the beam created is pulsed, not continuous, and the 
overall efficiency is low. 

Two new methods for creating narrow-velocity ground state atom barium beams 
have been developed in our laboratory. 

4.3.2 Atomic velocity selection via optical pumping 

We have developed a two-color optical pumping method for selectively preparing 
atoms in the ground state at a certain velocity, in an effusive beam. The beam is 
created in a cylindrical, externally heated oven made of Inconel. Continuous-wave 
lasers are used to pump tbe barium atoms, first into the metastable X D 2 state via 
the 6s 2 1 So *-> 6s6p x Pi cycling transition (A=553 nm), and then back to the 6s 2 1 So 
ground state with velocity selectivity via a tilted laser tuned to the 583 nm 6s5d X D 2 <-> 
5d6p 1 P\ transition (Fig. 14-1 ip . When the most probable velocity Vq is selected we 
obtain Av FW MM/t>o ~ 10% with approximately 50% repumping efficiency (Fig. 14-11]) . 
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Figure 4-11: Left: Twocolor optical pumping scheme. Right: Velocity profile of the 
selected atoms. In this experiment the pump laser is tuned to resonance, wherease 
the repump laser is tuned to resonance with the atoms with most probable velocity 
while the probe laser is scanned. In this example 5v(FW)/v ~ 13%. 



4.3.3 Supersonic oven 
Background 

Supersonic beam techniques have found many applications in spectroscopy and atom 
interferometry since their invention in 1951 [7H1 EH1 EH] • I n a supersonic beam, adi- 
abatic expansion as the beam exits the oven nozzle causes a dramatic cooling in the 
frame of the moving atoms or molecules; random thermal energy is converted into di- 
rected kinetic energy. In addition, the net flow rate is generally much higher than for 
an effusive beam. This solves the two problems associated with conventional effusive 
beams: low intensity and a wide distribution of velocities. 



Supersonic barium oven: apparatus 

We have developed a supersonic beam oven for barium closely following the design 
of Ref. [HI], which was based on an effusive beam source developed in our laboratory 
and previously used for the microlaser experiment [H]. 

The atomic beam oven consists of resistively heated, barium-filled piece of tan- 
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Figure 4-12: Relevant 138 Ba energy levels for optical velocity selection scheme. Dark 
arrows represent transitions used for pump (553.5 nm) and repump (582.6 nm). A = 
Einstein A coefficient; BF = branching fraction. 
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talum tubing, in which a small nozzle hole has been drilled. A diagram in shown in 
Fig. 14-131 The tubing is 6.5" long, 0.5" in diameter, with 0.015" wall thickness. The 
nozzle diameter varied from 150 jum to 406 /xm. 

Like the atoms strontium and ytterbium used in Ref. J^Tj , barium has two valence 
electrons and negligible dimer formation. 

The oven is created by the following procedure. A length of tubing is measured 
and cut with a pipe cutter. The bottom 1.125" is pinched flat in a vise between 
two aluminum blocks. One quarter of a 2 cm diameter barium rod (cut by bandsaw 
lengthwise) is inserted into the tantalum tube. The typical mass of barium in one 
oven is 12-16 g. The top 1.5" of the oven is then flattened in a vise, parallel with the 
bottom of the oven. The oven is mounted in a specially designed holder to position 
it with flat parts horizontal. It is clamped in a milling machine where a 1/8" end 
mill removes about 0.007" from the future location of the nozzle. This gives a shorter 
"channel" for the atom flow and creates a flat surface for drilling the nozzle, which is 
done in a miniature drill press. 

In the experiment, the oven is clamped by two screws at each end between copper 
electrodes. The upper electrode is connected to a copper tube which extends through 
an 8" Conflat flange via a 3/4" diameter Ultra- Torr fitting. The bottom electrode 
is soldered to a 4" diameter copper cold shield tube which both absorbs the oven's 
large thermal radiation and conducts current back up to the flange coaxially with 
the central conductor, to minimize stray magnetic fields from the oven. The upper 
and bottom electrodes and the cold shield are cooled by water flowing through 1/4" 
copper tubing which has been soldered to the copper pieces. 

A Proteus flowmeter placed on the return path of the oven cooling water is con- 
nected to a lamp which provides visual feedback that the water is flowing. 

The oven drive current is controlled manually using a variable Variac autotrans- 
former with 208 VAC line input and 0-230 VAC output. The Variac output is con- 
nected to a step-down transformer with 1:8 turn ratio, and the resulting output is 
connected to the center conductor and top oven flange (ground) via 2 pairs of 1" di- 
ameter welding cables. The oven voltage and current are monitored by AC voltmeter 
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Figure 4-13: Supersonic oven design 
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and a clamp-on ammeter. The voltage varies from 0-5VAC and current 0-1000 A 
rms?. The resistance of the entire circuit varies from approximately 4 mQ to 6 mfl 
depending on oven temperature. 

Supersonic barium oven: results 

As oven power is increased, the beam density increases and the beam undergoes 
a transition from Maxwell-Boltzmann (thermal) velocity distribution to supersonic 
velocity distribution. 

We measured the atomic density and velocity distribution via 553 nm 1 So <-> x Pi 
fluorescence. To measure density, the 553 nm beam was aligned to overlap the 791 
nm probe beam along the cavity axis and the fluorescence was detected by CCD 
(see section |4~6*|) . Velocity distribution was measured by aligning a focused 553 nm 
beam in the cavity in a direction opposite to atom beam direction, and detecting 
the Doppler-dependent fluorescence via a photon counting APD (EG&G C30902TE) 
previously aligned with the cavity transmission. A 500 fim diameter pinhole located 
in front of the APD and mounted on an XY translation stage blocks scattered light 
from the ground glass beveled edges of the cavity mirrors. 

Figure 14-141 shows the velocity distribution for the supersonic beam used in the 
experiments. This oven had a nozzle size of 406 /mi and was driven by about an rms 
current of about 600 amperes. 

The degree of velocity compression can be characterized by the ratio between 
FWHM of the measured velocity distribution and most probable velocity Av/vq. 

The measured velocity distribution f m (v) is equal to the actual velocity distri- 
bution f{y) convoluted with a lineshape function accounting for atom linewidth and 
isotope shifts. Laser power broadening effects can be minimized by using a probe 
laser beam intensity well below saturation. 

Although the supersonic expansion results in a higher velocity and thus smaller 
interaction time (roughly by a factor of 2), there are several advantages of the su- 
personic oven. First, and most importantly, velocity narrowing gives higher densities 
than velocity selection because non-selected atoms are lost in the latter. In addi- 
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Figure 4-14: Fluorescence signal of a counterpropagating 553 nm probe beam as mea- 
surement of supersonic beam velocity distribution. Zero of distribution determined 
by simulataneous measurement of Lamb dip transmission. Curve: Gaussian fit with 
peak velocity Vq = 816 m/s, a = 70.2 m/s. When natural linewidth and isotope shifts 
are taken into account, we find A^fwhm/^o = 16.5%. 

tion, the supersonic beam gives far higher densities than even a non-selected thermal 
beam, due to the higher oven temperatures. Experiments with the supersonic beam 
are simpler because the additional lasers and locking systems are not needed. Due to 
these advantages, all of the microlaser experiments in this thesis were done with the 
supersonic beam. 



4.3.4 Final aperture 

The small dimensions of the microlaser cavity mode place strong constraints on the 
size and location of the atomic beam. It must be small and closely centered on the 
mode if reasonably uniform coupling is to be obtained. 

Variation due to the Gaussian shape of the cavity mode can be made relatively 
small by limiting the dimensions of the atomic beam inside the cavity. For an atomic 
beam of height h centered in the mode, the ratio of minimum to maximum coupling 
is 

grain = e -(h/2w ) 2 (4.10) 

This value is plotted in Figure PFTol For a beam size of 25 //m we have g m i n / g m ax ~ 
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Figure 4-15: Atom beam width compared with cavity mode. Atom beam should 
"sample" only the peak of the mode to minimize coupling variation. Here the aperture 
width is set equal to one-half the mode waist w m . 



In order to create a small beam in the cavity we collimated the atomic beam 
with a 24.5 /im diameter aperture placed approximately 3 mm from the cavity axis. 
The distance between aperture and cavity is important because a larger distance 
increases the broadening of the atom beam due to the finite solid angle of the atom 
beam source (i.e. oven nozzle). The aperture was mounted onto a thin, hollow, 
tapered "pencil" which fit into the hole drilled into the PZT to allow passage of the 
atomic beam. The pencil is attached to the center of a 1/2" stainless-steel compact 
mirror mount which allows the orientation of the pencil to be modified precisely. 
The mirror mount was mounted on a miniature motorized rotation stage (National 
Aperture MM-3M-R) . A stepper driver (National Aperture MC-II) located outside 
the chamber is connected to the rotation stage by a 10-pin electrical feedthrough. The 
mirror mount adjustment screws were tweaked so that center of rotation is located 
at the center of the aperture. Finally, the rotation stage and aperture assembly were 
connected to an XY miniature translation stage (National Aperture ST1XY) which 
could be manipulated from outside the chamber through two 1/8 inch Ultra- Torr 



0.90 
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feedthroughs. The feedthroughs couple to the stage via flexible couplings. 

The pencil has a tapering design in order to be small enough at the tip to approach 
the cavity closely, and yet large enough to accomodate a range of atomic beam angles 
and avoid clipping the beam. The pencil was made from several pieces of brass 
telescoping square tubing, glued together with Varian Torr-Seal vacuum epoxy. The 
aperture (12.7 /im thick stainless steel, from National Aperture) was attached to the 
pencil tip with Torr-Seal. 

To center the atomic beam on the cavity mode axis, we monitored the microlaser 
emission power as a function of atom beam vertical position (measured by imaging 
the aperture with the imaging system described in Sec. 14.6(1 . By this method we 
are able to find the center of the mode to a few /zm . A typical result is shown in 
Fig. 14-191 Repeated measurements are taken as the aperture position is varied; low 
points are due to fluctuations in chamber pressure when the feedthrough is manipu- 
lated. Aperture position is measured by the CCD pixel number corresponding to the 
maximum intensity in an image of the oven blackbody radiation passing through the 
aperture. 
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Figure 4-17: Top view of aperture assembly including motorized rotary stage, mirror 
mount for angular adjustments, tapering "pencil", and final aperture in relation to 
optical cavity. Not shown: X-Y translation stages attached to rotation stage. 

4.3.5 Atom density modulator 

The density, peak velocity, and velocity compression of a supersonic beam are all 
functions of the oven current (i.e. temperature) and the amount of barium remaining 
in the oven. During our experiment, we would like to vary the density only. Therefore 
the oven current alone cannot be used to modulate density in the cavity, as it was for 
the original microlaser experiment [8] or some recent experiments jHl] 

We have introduced a variable beam density modulator which intercepts the 
atomic beam between the oven and cavity and prevents some or all of the atoms from 
entering the cavity. The modulator consists of a piece of sheet aluminum attached 
to a motorized linear feedthrough (Huntington L-2251-2) mounted at the 6- way cross 
between oven and cavity. This is shown in Fig. 14-201 The aluminum piece was ini- 
tially cut in a wedge shape in order to give a longer modulation range as it is moved 
vertically. This was later found to be unnecessary as sufficiently high resolution was 
obtained with a flat piece of aluminum. The modulator is located 21 cm from the 
oven; the cavity is 42 cm from the oven. 

The linear feedthrough has a lead of 0.05" per revolution. The stepper motor 
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Figure 4-19: Microlaser output for microlaser for constant density in initial threshold 
region, as a function of aperture position in vertical direction. 1 pixel corresponds to 
about 4.5 jum. 

turns one revolution per 200 steps, giving a linear resolution of 6.35/mi per step. 
Higher resolution is possible via microstepping. Full linear travel is 2" . 

The stepper motor is controlled by an Intelligent Motions Systems Panther LI2 
(aka Huntington MLC-11) Driver. This driver was found to create a great deal of 
electromagnetic interference in the laboratory, which was reduced using a set of opto- 
isolators and line filters. The driver is controlled by a PC via a standard RS232 serial 
connection. The driver is equipped with 3 logic-level auxiliary outputs which were 
used to control other experimental parameters (described in section l4.ll.lj) . 

The cavity fluorescence signal as a function of modulator position is shown in 



4.4 Laser systems, locking, and stabilization 
4.4.1 Laser equipment 

Two tunable lasers were used in the microlaser experiments. To provide the 791nm mi- 
crolaser pump beam and cavity probe, we used a Coherent 899-21 Titanium- Sapphire 
ring laser pumped by an Innova 310 argon ion laser operating at 8W in multiline vis- 
ible mode. 



Fig.E2H 
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Figure 4-20: Atomic beam density modulator 
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Figure 4-21: Beam density (fluorescence CCD counts per ms) as a function of beam 
modulator position; 1 step = 6.35 /zm. 

A Coherent 699-21 dye laser tuned to the 553 nm 1 So *-> x Pi transition was used 
to monitor atomic density in the cavity and measure the atomic velocity distribution. 
The dye used in this laser was Rhodamine 110. The laser was pumped by an 8W 
multiline visible beam from a Coherent Innova 200 Argon ion laser. 



4.4.2 Frequency stabilization of 791 nm beam by Doppler- 
free FM spectroscopy 

Frequency locking of the ThSapphire laser to the 1 So <-> 3 Pi transition was done by 
a standard FM modulation technique. This locking system is described in detail in 

An evacuated stainless-steel cylindrical cell containing about 2g barium is heated 
to 600° C by external heating elements to produce a barium vapor. Pump and probe 
laser beams pass through the cell from opposite directions, through two glass win- 
dows. The probe beam first passes through an electro-optic phase modulator (EOM) 
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Figure 4-22: Lamb dip used to lock Ti:Sapphire laser to 791 nm 1 So <-> 3 Pi transition. 

driven at 24.6 MHz. After passing through the cell the probe beam is incident on a 
photodiode (Thorlabs DET110). The signal is amplified and demodulated via mixing 
with the 24.6 MHz local oscillator signal in a double-balanced mixer (Minicircuits). 
The relative phase of the two inputs to the mixer is adjusted by varying the mod- 
ulation frequency and the lengths of BNC cables to produce a maximum dispersion 
signal, shown in 14-231 The dispersion signal is amplified and shifted to remove any 
DC offset from the mixer. It has a zero crossing at the peak of the Lamb dip (atom 
resonance). To lock the laser to the zero crossing, the dispersion signal is fed back to 
to Ti: Sapphire laser as an error signal in place of the signal from the reference cavity. 
The linewidth of the lock is estimated as approximately 300 kHz [83J. 

4.4.3 Locking of 553 nm beam to atomic beam 

The 699-21 dye laser was locked to the 1 So <-> x Pi transition via a different Lamb 
dip technique. After reflection through a beamsplitter, a strong 553 nm pump laser 
beam is incident on a barium cell heated to 490 C. The transmission through the cell 
is apertured then retroreflected back into the cell, overlapped by the original beam. 
The reflected beam forms the probe, which is detected by a photodiode (Thorlabs 
DET110) after the beamsplitter. The resulting absorption signal shows a Lamb dip 
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Figure 4-23: Dispersion signal used for locking of Ti: Sapphire laser to 1 So «-> 3 Pi 
transition at 791 nm. 



The laser locking system consists of an integrating amplifier and summing ampli- 
fier. A local oscillator at about 300 Hz frequency, 100 mV peak to peak amplitude 
dithers the frequency of the dye laser through the external horizontal input of the 
laser control box. The photodiode signal and local oscillator are connected to a lock- 
in amplifier. After adjusting the phase, the output signal is of dispersion type, with 
zero crossing at the peak of the Lamb dip. This signal is fed to an integrating ampli- 
fier, which controls the DC offset of the laser external input. The resulting lock stays 
within approximately 2 MHz of resonance. This is much smaller than the 19 MHz 
natural linewidth of the transition and was judged to be adequate. 

4.4.4 Intensity and position stabilization 

Both the 791nm and 553nm lasers experience some short-term and long-term intensity 
fluctuations. The dye laser has considerable ( ~ 20%) noise in the kHz range and the 
average intensity decreases gradually over hours as the dye decays. The Ti: Sapphire 
laser has little short-term noise but does have slow intensity fluctuations, most likely 
due to pump beam wander. 



(Fig.E23). 
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Figure 4-24: Lamb dip signal used for dye laser locking to the 553 nm 1 S <-> x Pi 
transition 
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Figure 4-25: 553 nm Lamb dip setup for 699 dye laser locking. 
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Figure 4-26: Laser intensity and position stabilization system for 791 nm and 553 nm 
beams 

Pump beam wander, thermal effects, daily power peaking, and other effects cause 
substantial variations in laser beam direction from day to day and within single days. 
Often these effects are corrected for by aligning a laser beam through two adjustable 
irises. However, since many aspects of the microlaser experiment require very precise 
beam positions (micron scale) and angles (milliradian scale), irises do not provide 
adequate repeatability. 

To reduce laser intensity noise and stabilize long-term power, and to stabilize laser 
positions, we have developed a laser stabilization system. Acousto-optic modulators 
(AOMs) are used to cancel intensity fluctuations. Coupling through single-mode 
optical fibers allows very consistent optical alignment of beams with respect to the 
setup. 

The intensity stabilization system uses a photodiode and integrating amplifier to 
provide feedback through an AOM VIDEO signal. A photodiode connected to a 
variable-termination, transconductance amplifier monitors the laser power after the 
single mode fiber as shown in Fig l4-26l The photodiode signal Vpd is compared with 
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Figure 4-27: Circuit diagram for intensity stabilization system 
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a target voltage Vr and the difference is passed to an integrating amplifier, the output 
of which is connected to the AOM VIDEO input (amplitude modulation). The output 
voltage range is trimmed by potentiometers to correspond to a range from minimum 
to maximum AOM transmission. 

Polarization maintaining (PM) single-mode fibers are used to stabilize the position 
of the 791 pump beam and the 553nm probe laser beams. We use lm long Newport 
SP-V fiber for the 791nm beam and a 5m Newport SP-? fiber for the 553nm beam. 
A schematic is shown in Fig. 14-261 

4.5 Pump beam 

The pump beam is supplied by the FM- locked Ti: Sapphire laser. The problem of how 
to pump the atoms correctly proved to be more complicated than expected, and led 
to the rediscovery of a method for inverting atoms by a defocused laser beam. 

The pump beam intersection with the atomic beam must be sufficiently small 
that upper state spontaneous decay does not significantly affect the final degree of 
inversion. For the same reason, the pump beam spot must be placed close to the 
cavity mode, but without overlap. 

Other critical parameters include the angle between the pump beam and atom 
beam direction (ideally exactly 90 degrees), total pump power, location of pump 
focus, uniformity of the pump beam intensity over the width of the atomic beam, and 
laser frequency jitter about the atom resonance. 

4.5.1 Optical setup 

As described in Sec.X the pump beam is coupled through a single-mode polarization- 
maintaining (PM) fiber for intensity and position stabilization. The emerging beam 
is collimated by an aspheric lens and passes through a beam steering apparatus to 
raise the beam height above the vacuum chamber. 

The pump beam steering apparatus is illustrated in Fig. 14-281 It consists of 
mirrors and a lens mounted on translation stages such that one stage controls the 
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Figure 4-28: Pump beam optical setup. Optics are fixed on two small rectangu- 
lar breadboards attached to separate translation stages to control pump focus and 
translation independently. 



pump beam focus and one controls the pump beam position (in the z direction), 
nearly independently. 

The lens is a single element, BK7 piano-cylindrical lens with focal length 400? mm. 
Tests using the CCD camera as a beam analyzer showed that a nearly diffraction- 
limited spot size of 35 fxm could be achieved at the focus. 

The pump beam is reflected down into the cavity by an 2" diameter, 45-degree- 
angled dichroic beamsplitter which has reflectivity approximately 90% at 791 nm and 
10% at 553 nm. A dichroic beamsplitter is used instead of a mirror in order to allow 
imaging of 553 nm fluorescence from above the beam splitter. In addition, the weak 
791 nm transmitted beam can be used to analyze the focus and beam quality of the 
pump. 

The window on the top of the chamber is an anti-reflection coated glass disc held 
in place by a custom metal flange and Viton O-rings. 

A coil-generated magnetic field of approximately 50 gauss parallel to the direction 
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of pump polarization (perpendicular to the cavity axis) splits the m-states by about 
50 MHz and ensures that only transitions between sublevels m = and m = occur. 

4.5.2 Pump beam parameters 

Here we describe the optimization of the pump beam parameters in more detail. 
Pump beam focus 

The pump focusing lens is placed approximately one focal length away from the cavity. 
To find the focus more precisely, we overlap the 791 nm pump with a resonant 553 nm 
probe beam. With the imaging system (Sec. I4.6j) in the top viewing configuration, we 
take pictures of the fluorescence of atoms in the probe field. By measuring the width 
of the fluorescence spot as a function of the pump focus translation stage micrometer, 
we can find the focus of the 553 nm beam. However, this is in general displaced from 
the 791 nm focus due to lens dispersion and different laser beam divergences of the 
beam before entering the pump beam setup. To address this problem we used the 
CCD camera, without lenses, as a laser beam analyzer for the 791 nm and 553 nm 
beams transmitted by the dichroic beamsplitter. The CCD position was positioned 
the same distance from the beamsplitter as the cavity is from the beamsplitter. The 
camera was moved by a translation stage to locate the positions of the 553 nm and 
791 nm beam foci. This gave the displacements of the two beam waists, typically 
about 1 cm apart. The pump focus was then adjusted by this amount to place the 
791 nm focus as close as possible to the atoms in the cavity. 

This procedure was not perfect: as described in Sec. 14.5.51 significant pump defo- 
cus was found to occur. 

Pump position 

Pump position can be monitored by a number of methods involving the 553 nm tran- 
sition and the imaging system. However, the most effective methods is to monitor the 
microlaser output power as a function of pump position. When the pump is placed 
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downstream of the cavity mode, no emission is seen. For pump beam overlapping the 
mode, a signal with unequal traveling-wave peak heights is usually observed. The rea- 
son for the unequal traveling-wave peaks is not clear. For pump beam upstream of the 
cavity mode, a signal is observed that is nearly flat, with a small slope corresponding 
to increasing atomic decay with pump-mode distance. 

In our experiments we aim for about 100 /im separation between the pump field 
and cavity mode. At this distance, which is admittedly somewhat arbitrary, atom 
decay between pump and mode is quite small, < 10%. 



Pump angle 

To minimize Doppler broadening it is important for the pump beam to be as close as 
possible to perpendicular to the atomic beam direction. Two methods for ensuring 
this condition were explored. First, we overlapped a 553 nm probe beam over the 791 
nm pump as closely as possible, and blocked the 791 nm beam. Then we measured 
the linewidth of the 791 nm fluorescence as the dye laser was scanned in frequency, 
and repeated this process as the angle between laser and atom beam was varied. 
The linewidth increased linearly with any angular deviation from a right angle. This 
technique was very time-consuming and gave a pump beam perpendicular to the atom 
beam only to the accuracy that the 791 nm and 553 nm beams were parallel in the 
cavity. The second technique was to monitor the microlaser output for a constant 
atom density a few times over threshold, while the frequency of the 791 nm pump laser 
was varied over a 50 MHz range via the external input of the Ti: Sapphire laser. In this 
manner, it could be determined if the microlaser power was maximized for positive 
or negative detunings relative to the atom frequency (i.e. the atoms experienced a 
blueshifted or redshifted pump field). The pump angle was adjusted until the the 
microlaser power was maximum with the laser on exact resonance. 
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4.5.3 Shelving experiment 

In order to test the efficiency of the pump process, we perform a "shelving" experiment 
as follows. 

A 553 nm probe beam is aligned with the 791 nm probe which had been matched 
with the TEM 00 mode of the cavity. Therefore the 553 nm probe field is close to 
centered on the cavity mode (the beam size may be different, but it will be close to 
the focus since the 553 nm probe passes through the same achromatic lens). Using 
the CCD camera, we measure the fluorescence from the atoms in the cavity with and 
without the pump beam. Since the fluorescence is a measure of the ground state 
population, the ratio between the fluorescence intensities with and without the pump 
beam is equal to the average upper state population. If the pump beam were perfectly 
efficient, no fluorescence at 553 nm would be seen, because all atoms are in the upper 
state of the 1 S <-> 3 Pi transition. 

Shelving ratios as low as 10% (i.e. 90% of atoms in the upper state) have been 
measured. 



4.5.4 Pump power dependence of inversion 

When the first shelving data was collected, a puzzle presented itself. Atoms interact 
with the pump beam for a certain time dependiing on their velocity v, 

U = 0FWpUmP (4-11) 

v 

where w pump is the waist of the pump beam, assumed to be Gaussian. The excited- 
state probability for a given atom after leaving the pump beam is expected to be 

p(l) = gin 2 ^pump ^j (4 

Averaging over atom velocity distribution f(y) gives 

P. = jf pe£fi=) /(,) *, (4.13) 
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Figure 4-29: Ground state amplitude of atoms, measured by 553 nm fluorescence: 
monotonic dependence on pump power. Crosses (x): with B field. Circles(o): with- 
out B field. 

We expect to observe an oscillatory excited state probability damped by the ve- 
locity distributions. In the limit of strong field the excited state population should 
be 1/2. 

Oscillations as a function of pump power were initially not observed. Instead, we 
observed the dependence drawn in Fig. 14-291 The excited state population increases 
monotonically to approximately 80% and does not approach 1/2. 

4.5.5 Adiabatic inversion by defocused pump beam 

The unusual behavior was explained by the following effect. Suppose there is some 
distance y between the atoms and the focus (minimum waist) of the pump beam - 
for example, the beam comes to a focus before reaching the atoms. Then the atoms 
pass through a laser beam which is diverging in the direction of their propagation. 
As a result, the angle between atom velocity v and the wave vector of the pump field 
k varies approximately linearly with position, causing a Doppler-induced frequency 
chirp through resonance. This leads to an adiabatic inversion process as the atom's 
Bloch vector precesses about a relatively slowly varying "torque" which tends to align 
the atomic Bloch vectors with the inversion axis. A quantitative analysis follows. 
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Without loss of generality we restrict our attention to the x — y plane containing 
both the atomic beam and Gaussian beam axis. The phase function of a Gaussian 
beam depicted in Fig. 14-301 is (see e.g. 



y) = -ky + taa _1 (y/j/ ) - (4.14) 

where k is the wave vector, yo = ttw^/X is the Rayleigh range and the wavefront 
radius of curvature 

R(y) = y (i + f 2 ) (4.15) 

The angle 9 between the Gaussian beam axis and the normal to the phase front 
can be found via 9 = tan -1 (dy/dx) where dy/dx is the slope of the wavefront (i.e. 
line of constant $): 



9 = tan" 1 



,2\ i yo i x 2 y 2 -y\ 



-x 



-(y 2 + y 2 ) + f + 



k 2 y 2 +yf ) 



(4.16) 



The sign of the angle is chosen to be positive when the x component of the vector 
normal to the phase front is opposite the direction of the atom velocity. For simplicity, 
we consider only the case of atomic velocity perpendicular to Gaussian beam axis. 

The equation for 9 simplifies assuming k ^> y and x<i/, both valid easily in our 
case: 

9 « —^- 2 x (4.17) 
The Gaussian waist of the beam as a function of y is given by 

w(y) = w py f\ + (y/y ) 2 (4.18) 

An atom traversing the pump beam at vertical position y experiences, as a function 
of position, the following Doppler shift: 

I y\ 

A d (x) = kv sin# rs kv — ^x (4.19) 

V 2 + Vo 
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Figure 4-30: Atoms displaced from pump beam focus experience frequency chirp and 
adiabatic inversion. 
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The excited state may decay back to the ground state or to a metastable D state. 
The decay rate to the ground state is T g = 50kHz and to the metastable state is 
T m = 70kHz. The optical Bloch equations describing the pumping process including 
loss are 



Pee® = -^-R 2 (t) -T g + T m ) Pee (t) 



R 2 (t) = A(t)i?i(£) + n (t)( Pee (t) - Pgg {t)) - -fR 2 (t) 
R^t) = -A(t)R 2 (t) - ^-Ri(t) 



(4.20) 

(4.21) 
(4.22) 
(4.23) 

(4.24) 



where p ee , p gg , and p mm are excited state, ground state, and metastable state proba- 
bilities; 

R^t) =p eg + Pge (4.25) 



R 2 (t) = i(p ge - p eg ) 



(4.26) 



and A is the atom-laser detuning 



A = uji + A d - uj a 



(4.27) 



where io\ and u a are laser frequency and atom resonance frequency. The on-resonant 
Rabi frequency is given by 



Q (t) = n exp[-(vt/w(y)) 2 ] = J —T g exp[-(vt/w(y)) 



(4.28) 



where Qq and / are the Rabi frequency and on-axis laser intensity and I s is the 
saturation intensity 

■nhr 

(4.29) 



I* = ^-T g w 13/iW/cm 2 



s 3A 3 

Fig. 14-311 shows the results of numerical integration of the Bloch equations for different 
displacements y. For y = 3y the Rabi oscillations are almost completely absent, and 
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inversion occurs almost independently of pump intensity for / > 27. 
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Figure 4-32: Ground state probability as a function of ^Jp/p-k/2, where p is the pump 
power and p n /2 ~ 1.8/xW. 

If adiabatic inversion is responsible for the pump power dependence, we should be 
able to observe Rabi oscillations by focusing the pump beam closer to the atoms (i.e. 
decrease y). The results are shown in Fig. 14-321 Note that peaks and valleys of the 
oscillations occur at integer values of ^p/Pn/2- The curve is not asymptotic to 1/2 
for large pump fields; due to uncertainties in the experiment this data may reflect a 
intermediate state between pure Rabi oscillation and adiabatic inversion. 

A quantitative comparison between experiment and theory will be performed in a 
future study. It has come to our attention that adiabatic following due to a diverging 
laser beam was previously observed in [86] and |87j . These experiments involve atoms 
with more complex level structures than the present one. 
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Figure 4-33: Schematic of imaging system 



4.6 Imaging system 



A cavity imaging microscope was developed to perform several functions, described 
below and in following sections. 



4.6.1 Setup 

The optical design of the imaging system is shown in Fig. 14-331 It can be thought of as 
a relay lens system combined with a CCD microscope. A 2:1 lens system consisting of 
two 1" achromatic doublets creates a real image of the cavity. A microscope objective 
relays the image to a CCD camera. The CCD is positioned at the field focus plane 
of the objective. The use of microscope objectives gives well-corrected imaging and 
allows the magnification of the system to be varied easily without realignment. In 
our experiments we use the 5x objective for the side imaging and lOx for the top 
imaging (see below for descriptions of the two imaging configurations); the lOx and 
20x objectives can also be used in the side imaging configuration when the position 
of the final collimation aperture needs to be determined precisely. 
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Figure 4-34: Side and top viewing imaging system configurations. Lab jacks and 
supporting structures not shown. 

The achromatic lenses are translatable in 3 directions by translation stages for 
centering and focusing of the image. The optics, translators, and camera are mounted 
on an aluminum plate which is fixed to a laboratory jack. 

For the side imaging configuration, a 1 inch diameter interference filter (Oriel 
53890, 550 ± 5 nm) is placed between the two achromats to block transmission of the 
oven's blackbody radiation. For the top imaging configuration, a color glass filter is 
used to eliminate pump beam reflection and scatter into the imaging system. This 
filter has transmission T > 95% at 553 nm and T <0.1% at 791 nm. 

Background light is blocked by a plastic C-mount coupled tube connecting the 
CCD camera and objective lens, and a sliding cardboard tube covering the optical 
path between the objective and relay lens. 

The CCD used is a Roper Scientific / Photometries Sensys 400. This camera uses 
a thermoelectrically cooled, 1/2-inch format, 768 x 512 pixel Kodak KAF 0400 CCD. 
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configuration 


lens 1 f.l. 


lens 2 f.l. 


objective (typ.) 


magnification 


side view 


20 cm 


10 cm 


5x 


2.5 


top view 


45 cm 


10 cm 


lOx 


2.2 



Table 4.1: Optics for two imaging system configurations 



Readout digitization is 12 bits. 

The camera can be used in two configurations (Fig. l4-H4l In the first configuration, 
the camera views the cavity from the side, antiparallel to the atomic beam direction. 
In this view the aperture can be imaged and fluorescence of the atom beam appears as 
a small spot. In the second configuration the camera images the cavity from above. A 
mirror is mounted at a 45-degree angle above the top window of the vacuum chamber, 
bending the optical path so that the optics and CCD can be mounted horizontally, 
above the optical path for the side-viewing configuration. The optics used for the two 
configurations is summarized in Table l4~Tl The CCD mount is designed so that the 
CCD can be swapped between the two configurations easily. 

4.6.2 Uses of imaging system 
Atom-cavity alignment 

To achieve the most uniform atom-cavity coupling possible, it is essential that the 
atomic beam be centered on the cavity mode. This can be done by adjusting the 
final aperture vertical translation stage to maximize output power. However, the 
miniature translation stage is manipulated via a flexible coupling through an Ultra- 
torr vacuum fitting and exhibits considerable backlash (hysteresis). In addition, the 
microlaser output power is not a highly sensitive measure of the centering of the 
aperture on the mode, especially when the photon number stabilization effect is taken 
into account. Therefore is it useful to be able to monitor the aperture position 
accurately. This can be done by focusing the imaging system on the oven's blackbody 
radiation emitted from the aperture. A computer program monitors the microlaser 
output and blackbody peak vertical position as the latter is varied manually. 
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Figure 4-35: Atom fluorescence imaged by microscope in side imaging configuration. 
Similar but smaller images were used to measure atomic density during the experi- 
ment. 

Density calibration 

The density of atoms in the atomic beam, and therefore the average number of atoms 
in the cavity, can be monitored quantitatively by detecting fluorescence of atoms in 
the cavity. Details of this experiment are given in Sec. 14.91 

In addition to density calibration, the imaging system can be used to monitor the 
condition of the final atom beam aperture. Over time this 25 /zm aperture becomes 
coated and eventually clogged with barium. By monitoring the fluorescence and 
blackbody radiation, we can tell when the aperture needs to be removed and cleaned. 

Characterization of laser and atomic beam 

In order to calibrate density, focus the pump and probe beams, and position the 
pump beam correctly, it is necessary to determine the position and focus properties 
of the 791nm and 553nm laser beams. The imaging system can be used to measure 
the focus and position of these beams. 

Shelving measurement 

The imaging system was used to detect 553 nm fluorescence in the cavity while the 
791 nm pump beam was present. This allows us to measure the degree of excitation 
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Figure 4-36: Microscope image of center of cavity taken from side imaging configu- 
ration. Left and right: cavity mirrors. Center: final collimating aperture. For this 
image, cavity was illuminated by flashlight and oven blackbody radiation was used 
to image aperture slit. 

due to the pump beam. 

4.7 Detectors 

4.7.1 Avalanche photodiodes 

A photon-counting, thermoelectrically-cooled avalanche photodiode (EG&G 10902TC) 
is used to detect the 553 nm fluorescence in the cavity when measuring beam velocity 
distributions. An uncooled APD (EG&G f0902) detects the cavity transmission for 
cavity locking and finesse measurements. 

4.7.2 Photomultiplier tubes 

A photomultiplier tube was used to detect microlaser emission. Our PMT is a Hama- 
matsu R943-2. It is cooled to less than —20° C by a thermoelectric cooler (Pacific 
Instruments 3470). At this temperature it has a dark count rate of about 20 cps. The 
bias voltage used was -1 500V. The signal from the PMT is terminated into 50 f2 and 
discriminated/counted by a Stanford Research SR400 photon counter. 
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Figure 4-37: Photo of the microlaser experiment. Atom beam oven is attached to top 
flange of vertical cylinder portion of chamber at right. Chamber containing cavity is 
at center. 

The PMT count rate was verified to be linear for the count rate used in the 
experiment. The saturation count rate was not measured and may be significantly 
higher than the damage threshold of the PMT. 

4.8 Vacuum system 

A vacuum chamber houses the cavity and atomic beam oven. It consists of two large 
cylindrical sections, the first containing the oven and the second containing the cavity. 
The two chambers are connected by flexible welded bellows couplings. The atomic 
beam passes through the flexible segment and into a 6-way cross before entering the 
cavity chamber. A diagram is shown in Fig J4-38l 

All ports are faced with Conflat-type flanges. OFHC copper and Viton gaskets 
are used to create seals. The oven chamber is connected to a diffusion pump (Varian 
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beam blocker 

Figure 4-38: Diagram of major chamber components showing atom beam path. 

VHS-4) with pumping speed 1000 1/s. A pneumatically driven 8" gate valve (MDC 
GV-50000) separates the oven from the chamber. An interlock circuit ensures that 
the diffusion pump can only be operated when the chamber and foreline pressures 
(measured by two thermocouple gauges) are sufficiently low, and the diffusion pump 
cooling water is flowing. The foreline of the diffusion pump is connected by hoses to 
a roughing pump located in an adjoining room to minimize vibrations and noise. The 
foreline pressure was about 5 millitorr, measured by a thermocouple vacuum gauge. 
A tubulated ion gauge measures residual pressure in the chamber, which during the 
experiment is typically ~ 10~ 6 Torr. 

The vacuum system was fabricated by MDC Vacuum Products, and several mod- 
ifications were performed by Sharon Vacuum and the MIT Central Machine Shop. 

From time to time it was necessary to bring to chamber to atmospheric pressure, 
in order to change the oven, reload the barium cell, etc. In these cases we filled the 
chamber with argon to avoid oxidizing the barium in the chamber. This also avoids 
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Figure 4-39: Gaussian mode and equivalent tophat profile 



introducing water vapor into the chamber which is difficult to pump out. 

4.9 Density calibration 

We now describe the experiments and calculations used to determine the number 
of atoms present in the cavity. First, let us clarify the definition of the intracavity 
atom number N. Since the cavity mode is Gaussian, it is not clear where a boundary 
between the inside and outside of the mode can be drawn. We define N to be the 
average number of atoms within the boundary of a top-hat distribution (Fig. I4-39|) 
with height equal to that of the Gaussian and with equal area (i.e. total interaction). 
This condition is met for a top hat of width w top hat = \f^w m . 

Note that we assume that each atom in the atom beam passes through the mode at 
some time; this is true because the atomic beam size in the direction perpendicular to 
the cavity axis is smaller than the mode waist. By contrast, in the original microlaser 
experiment the atomic beam was 300 fim in diameter, much larger than the mode. 

To measure the atomic density we overlap a 553 nm probe beam onto a 791 nm 
beam which has been focused to match as close as possible only the TEM 0Q mode 
of the cavity. The 553 nm beam then intersects the atomic beam perpendicularly 
except for any cavity tilt that is introduced. The resulting 1 S <-> x Pi fluorescence is 
collected by the camera in the side viewing configuration. 
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Figure 4-40: density- measurement-setup 



The total number of fluorescence photons incident on the CCD is given by 

n p = VopticsFnW F s T ex (4.30) 

where Coptics is the transmission of the optical system including chamber window, 
Fn is the fraction of scattered light which travels in the direction of the collection 
optics, and W is the total scattering rate. T ex is the exposure time. Fg < 1 is a 
detuning factor reflecting the Doppler broadening due to the tilted cavity, and is 
approximately equal to the ratio between the homogeneous and inhomogeneous atom 
linewidths. 

The transmission Coptics is the product of the transmission of the chamber window, 
two lenses, 553 nm interference filter, and microscope objective. This was measured 
by passing a laser beam through these optics: 

?? optics « 22% (4.31) 

The parameter Fq involves the solid angle of collection and the dipole radiation 
pattern of the atoms. The polarization of the probe field is set to be vertical so that 
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the camera collects light scattered at the optical 90 degrees from the dipole direction. 
The fraction of power emitted into the solid angle subtended by the camera optics is 



/optics Sm (6)dQ ^ ^optics _ 3 Coptic 

J all sin 2 (6)dtt ~ 2njfd6 sin 3 (6) ~ 2 4vr 



(4.32) 



where we have assumed 6 does not deviate significantly from tt/2 over the solid angle 
of the collection optics. This result simply means that the fluorescence detected at 
90 degrees from the dipole direction is 50% than that expected from an isotropic 
radiation pattern. 

The solid angle of the optics is that of a 1" lens located 20 cm from the fluorescence 
source 

tt(0.5) (2.54cm) 2 , A . 

Optics « 1 ( 2Q cm)2 = °- 019 ^ ( 4 - 33 ) 

The total rate of photon scattering from the source is 

W = J drp(r)R(r) (4.34) 

where p(r) is the atom density and R(r) is the scattering rate per atom at point r: 

P(r) = Fa / ( :r )/ /sat r a {p/p s )exp{-2{x/w p ) 2 ) 
[) 2 1+/(i)// 8at 2 l + (p/ Ps )eM-2(x/w p )i) { - 1 

Here T a m 19MHz is the natural linewidth of the 1 So <-> : Pi transition; w p is the 
mode waist of the probe beam; I s is the saturation intensity; p s is the power of the 
553 nm probe, measured outside the chamber, such that the laser intensity center of 
the focused probe beam is I s . We then have 

W = T dx\'^ {p,Ps) ex P(~ 2 (*M) 2 ) (4 36) 

J-oo 2 1 + (p/p s ) exp(-2(x /w p ) 2 ) v ' 1 

where A' = / dy dz p is the constant linear density along the atom beam direction. 

The parameter p s is measured by using the camera in the top viewing configu- 
ration and taking exposures for different probe powers. The peak of the measured 
fluorescence distribution corresponds to the most intense (axial) part of the probe 
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Figure 4-41: Saturation curve of peak of 553 nm fluorescence: fluorescence counts vs. 
input probe power, fiW. 



beam; by plotting this against the probe power and fitting a saturation curve we 
obtain a saturation power (Fig. 14-4 lj) of 



p s = 11.2/xW 



(4.37) 



The probe beam waist is measured with the same camera data by imaging the 
fluorescence distribution directly (Fig. 14-42)) . After accounting for atom saturation 
and camera defocus/aberration effects we obtain 



w p = 139/im. 



(4.38) 



The total integrated counts recorded on the CCD is given by 



C = n p ■ r] CCD /G « 0.0625n p 



(4.39) 



where r/ccD ~ 0.35 is the CCD quantum efficiency at 553 nm according to manu- 
facturer specifications, and the gain is G ~ 5.6 electrons per ADU (analog-to-digital 
unit) at the camera's gain index 3 setting. 

110 



fluor. cts 



15 



12 . 5 



2.5 



7 . 5 



10 



5 



v_ 



*25 




100 12-5-' 150 



50 



75 



Pixel number 

Figure 4-42: 553 nm fluorescence distribution from CCD. Fit includes saturation 
effects and optical aberration function. 

The CCD exposure time was typically set in software to be 300 msec. However, it 
was discovered (Fig. I4-43J) that the total fluorescence measured was not proportional 
to the exposure time when those exposure times were small; this is probably due to 
finite shutter opening and closing times. An extra factor was used to compensates 
for this effect: T ex = 142 ms for a set time of 300 ms. 

The lineshape of the 553 nm transition was measured by varying the laser detuning 
from the computer and recording the total fluorescence. The resulting curve was fit 
accurately by a Lorentzian with FWHM of 34 MHz. The detuning correction factor 
is then F 5 = 19MHz/34MHz = 0.559 

Finally we calculate the intracavity atom number: 



The ratio between CCD counts and atom number is most conveniently expressed 
as the number of counts per millisecond of exposure time, per atom, with probe power 
10 fiW. We find that N = 1 atom corresponds to 1.16 counts per ms on the CCD; 1 
cpms corresponds to 0.862 intracavity atoms. 



(4.40) 
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Figure 4-43: CCD counts per millisecond (cpms) for a constant atom density when 
the exposure time T is varied. Drop-off for short exposures is probably due to opening 
and closing times of mechanical shutter. 



4.10 Photon number calibration 



The PMT count rate n d as a function of cavity photon number n is 

T 

n d = n( + A + S) ) r ^°ptics^AOMr/PMT (4.41) 

where T = 0.5 x 10~ 6 , A = 0.2 x 10" 6 , and S = 3.0 x 10~ 6 are estimates of the 
transmission, absorption, and scattering coefficients of the cavity mirrors; the factor 
of 2 accounts for the collection of light from only one side of cavity. T c is the cavity 
loss rate; ^aom is the fraction of light switched by the AOM after the cavity into 
the PMT (i.e. diffraction efficiency, about 70%); Coptics is the transmission of optics 
including the rear mirror surface, chamber window, several mirrors and lenses guiding 
the beam to the PMT, and the PMT window. ?7pmt is the quantum efficiency of the 
PMT at 791 nm, about 12% according to manufacturer specifications. 

Our estimate for the number of photons in the cavity is n = 0.16 per 10 3 counts 
per second on the PMT, with no filter in front of the PMT. 

Both the estimate for atom and photon number are subject to fairly large system- 
atic error; each may have errors as large as a factor of 2. 
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4.11 Data acquisition and computer control 



The experiment is controlled by a Gateway 200 MHz PC running Windows 95. The 
sequencing and data collection routines are written in VPascal, the language used by 
the camera control software, Digital Optics V++. The sequencing programs used in 
this experiment can be found in Appendix IU1 

4.11.1 Experiment automation 

Table I4~21 summarizes the computer automation of various devices. These are further 
described and commented on below. 



Device 


Level 


Control 


SR400 photon counter 
IMS atom beam mod. driver 
Laser intensity mod. (791 nm pump) 
Laser intensity mod. (553 nm probe) 
Analog multiplexer (MUX) 
Cavity lock disable or scan trigger 
791nm probe beam AOM Vr 
AOM after cavity V T 


digital data 
digital data 
0-1 V 
0-1 V 
TTL x 3 
TTL 
10V(on), 4V(off) 
10V(on), 4V(off) 


PC serial port 2 
PC serial port 2 
SR400 PORT1 output 
SR400 PORT2 output 
IMS TTL OUT 1-3 
MUX output 1 
MUX output 2 
MUX output 3 



Table 4.2: Summary of computer automation setup 



Laser intensity modulation 

The two analog outputs of the SR400 photon counter are used to control the intensities 
of the 553 nm probe and 791 nm pump. During the "off" states the target voltage is 
set to -0.1 volt in order to ensure that the minimum intensity is reached. 

Atom density modulation 

The stepper motor driving the density modulator is controlled by the computer 
through the serial port. The stepper rate was 400 steps per second. No microstepping 
was necessary. 
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AOM switching 

Two acousto-optic modulators required binary control: the 791 nm probe switching 
between cavity locking and data collection, and the AOM deflecting the cavity emis- 
sion to the setup during data collection. Both AOMs were controlled by TTL 
signals from the stepper driver. 

4.12 Cavity PZT control 

The cavity frequency may be scanned through the microlaser resonance or locked onto 
it. In some sense a scanning experiment is more general since it includes variation in 
atom-cavity detuning. Other experiments require the laser to be kept on resonance 
while other parameters are varied. 

Cavity scanning 

For cavity scanning, the PZT is simply driven by a ph 1 V ramp signal generated 
from an oscilloscope, which sweeps in AUTO' mode or is triggered by a pulse from 
the TTL output from the computer. To create bidirectional scans we designed an 
additional circuit which outputs |V; n (t) — Vol given an input voltage V in (t) from the 
oscilloscope or other ramp source. Vo is equal to half the peak-peak height of the 
ramp signal. The scan and sample PMT signal is shown in Fig. 15-11 

The experimental sequence for cavity scanning experiments is the following: 

1. Move atom density modulator for the desired density 

2. Turn off the 791 nm pump beam, 

3. Turn on the 553 nm probe, with power 10 /iW 

4. Measure atomic density by exposing the CCD in the side viewing configuration 
for a set time, usually 300 msec, 

5. Turn off 553 nm probe and turn on 791 nm pump beam 
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Figure 4-44: Probe transmission spectrum for frequency calibration of PZT scans. 
Sideband frequency shift from carrier is 24.67 MHz. 

6. Trigger oscilloscope to drive cavity PZT with a bidirectional ramp of duration 
about 4 seconds 

7. Start photon counter to collect data from PMT during cavity scan 

8. Download PMT count data to computer via serial connection 

It was occasionally necessary to adjust the HV PZT offset voltage manually to 
compensate for a slow PZT drift which would otherwise cause the microlaser reso- 
nances to drift out of the scan range. These adjustments were made between data 
collection intervals so that the linearity of the PZT scans was preserved. 

To calibrate the frequency of the scans we measured the cavity transmission of a 
very low power probe beam with no atoms present; the result is shown in Fig. 14-441 

The background counts level, due primarily to pump scattering from cavity mirrors 
and stray room light, is subtracted by computer. 

Cavity locking 

To maintain the cavity on resonance, we use the following cavity lock technique. 

A probe beam from the Ti:Sapphire laser is upshifted 100 MHz and then down- 
shifted about 113 MHz by a pair of AOMs. The frequency of the second AOM 
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Figure 4-45: Cavity locking circuit 



adjusted via the input tuning voltage Vr to match the peak of one of the microlaser 
resonances, in our case the lower frequency resonance at u a — kv^O. The AOMs are 
also used to switch the probe beam on and off. 

For cavity locking, a probe beam of about 1 /xW is incident on the resonator. The 
cavity transmission is focused onto an APD. A circuit compares the APD voltage 
with a fixed reference voltage adjusted to a value approximately one-third of the 
peak signal height. The voltage difference is integrated and the result fed back to the 
PZT control voltage. If no strong mechanical or acoustical disturbances occur, this 
system can remain locked for more than 1 hour. 

Experiments which involve cavity locking require us to switch between locking 
and data collection. The cavity can be unlocked for a maximum of about 1 second 
before the cavity drifts too far from resonance for lock to be regained. During cavity 
non-lock, a TTL signal applied to the lock box causes the integrator to hold its value 
by grounding its input. 
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Chapter 5 

Results and Analysis 



We describe the results of microlaser experiments with a high-density supersonic beam 
and nearly uniform atom-cavity coupling. 

Experiments were performed with the goal of answering the following questions: 
(i) Does the theory of Chapter |21 apply to the microlaser with N c q ~ 1000? (ii) Do 
multiple thresholds occur, and if so, where? (iii) What are the time scales of the 
transitions, and can hysteresis be observed? 

5.1 Study of microlaser with variable detuning 

We first measured the microlaser output as a function of atom-cavity detuning by 
applying a bidirectional ramp signal to the cavity PZT. The experiment sequence is 
given in Section fa. 121 The atom density modulator position was varied between 1350 
and 1730 with a step size of 10, and 1450 to 1500 with a step size of 2 (c.f. Fig |4-20|) . 

5.2 Cavity scanning data 

Data for the detuning curves for a range of atom densities are shown in Figs. 15-51 to 
15-241 The PMT data is folded back onto itself to show both results from the PZT 
scan in both directions. The solid and dashed lines represent scans in the positive and 
negative detuning directions, respectively. PMT counts are measured in time bins of 
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Figure 5-1: Bidirectional PZT voltage scan shown with typical PMT data for a high 
atom density. Horizontal range corresponds to 4 seconds. 



0.01 seconds. 

Effective atom numbers are following the estimates given later in this chapter. 
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Figure 5-4: Fluorescence 80.2 cpms (N eS = 17.6) 
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Figure 5-5: Fluorescence 100 cpms (N eS = 22) 
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Figure 5-6: Fluorescence 659.1 cpms (N eS = 145) 
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Figure 5-7: Fluorescence 950 cpms (N eS = 209) 
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Figure 5-8: Fluorescence 1850 cpms (N eS = 407) 
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Figure 5-9: Fluorescence 2393 cpms (iV cff = 526) 
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Figure 5-10: Fluorescence 2748 cpms (iV eff = 605) 
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Figure 5-11: Fluorescence 2889 cpms (N cS = 636). Second threshold is observed. 
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Figure 5-12: Fluorescence 2986 cpms (N cS = 657) 
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Figure 5-13: Fluorescence 3014 cpms (iV eff = 663) 
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Figure 5-14: Fluorescence 3075 cpms (iV e ff = 677) 
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Figure 5-15: Fluorescence 3252 cpms (iV eff = 715). Spikes are the first signs of third 
threshold. 
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Figure 5-16: Fluorescence 3432 cpms (N eS = 755) 
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Figure 5-17: Fluorescence 3556 cpms (iV eff = 782) 
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Figure 5-18: Fluorescence 3595 cpms (iV eff = 791) 
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Figure 5-19: Fluorescence 3628 cpms (N eS = 798) 
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Figure 5-20: Fluorescence 3690 cpms (iV efr = 812) 
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Figure 5-21: Fluorescence 3996 cpms (N eS = 879) 
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Figure 5-22: Fluorescence 4008 cpms (N cS = 882) 
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Figure 5-23: Fluorescence 4585 cpms (iV efr = 1009) 
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Figure 5-24: Fluorescence 5140 cpms (N cS = 1131) 
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As density is increased from zero, the two resonances broaden and increase in 
amplitude. At about 3000 cpms (Fig. I5-11JI a "spike" appears near the peak of each 
resonance. We will show below that this is a transition to the second branch. With an 
further increase in density to 3250 cpms (Fig. 15-15(1 a jump to another, third branch 
appears. 

Above second threshold, the positive and negative direction scans are in general 
quite different in shape, although they overlap over some intervals. Transitions in the 
positive-detuning direction occur at equal or greater (more positive) detunings than 
the corresponding transitions in the negative-detuning direction. This is as expected 
from hysteresis behavior. 

To measure the frequency splitting between the two traveling-wave peaks in the 
microlaser, we found the frequency shift between the two maxima for each scan below 
the second threshold; the average value was 27.1 MHz. The measured peak of the 
velocity distribution was Vq = 815 m/s; therefore we infer the cavity tilt angle is 
9 ~ 13 mrad. 

5.2.1 Photon number vs. atom number 

The detuning curve data was analyzed by computer to give, for a fixed detuning, the 
microlaser emission rate as a function of atomic density. 

Data for center of the peaks A = ±13.5MHz is given in Fig. 15-251 in terms of 
an estimated cavity photon number n and effective atom number N e s. The curve 
is the result of the rate equation model of Chapter |2] with modifications for velocity 
distribution, nonuniform coupling, and detuning effects. 

Calibrations for the atom density and photon number were adjusted to give the 
best fit between theory and experiment as follows. 

The photon number calibration from Chapter 0] was n = 0.16 for every 10 3 cps 
on the PMT, with no filters. With the OD 1.0 filter used in this experiment, this 
becomes n = 1.6 for every 10 3 cps on the PMT. The experimental best fit is with 
n = 0.25 x PMT counts per 0.01 second bin, or n = 2.5 for every 10 3 cps. The 
experimentally derived calibration is 56% higher than the original estimate. 
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The estimate of atom number from Chapter 0] was iV = 1.16 for every count 
per millisecond (cpms) on the CCD. For the 300 msec exposures in this experiment, 
the measured fluorescence count rate is underestimated by a factor of 2.11 due to 
the nonuniform exposure effect (Fig. 14-43*]) . Therefore the expected value is N = 0.55 
atoms per cpms. The effective atom number is N e g = (0.80 — 0.20)iV = 0.33 per cpms. 
The experimental fit gives N e g = 0.22 per cpms, 33% lower than the estimate. This 
fitting procedure is justified by the large uncertainties in our estimates for absolute 
photon and atom numbers. 

Note that in the detuning curves, the positive-direction data and negative-direction 
data are in most cases nearly mirror images of one another. This is expected from 
the symmetry between the two traveling-wave components. We therefore choose to 
display the photon number vs. atom number data in following way: Circles (o) rep- 
resent data from (i) positive-direction detuning scans with positive detunings and 
(ii) negative-direction scans with negative detunings; crosses ( x ) give data from (iii) 
positive-direction detuning scans with negative detunings and (iv) negative-direction 
scans with positive detunings. 

The agreement between theory and experiment is, overall, quite good. A fairly 
large number of "outliers" are present. These points represent the system undergoing 
transitions between the first and second branches and reflect the fact that the second 
threshold occurs very close in detuning to the central point. To support this point, 
Figs. 15-261 and 15-281 show the microlaser output curves for nearby detunings 12 MHz 
and 15 MHz. Dramatic changes occur for very small changes in detuning. The line is 
the same as in Fig. 15-251 In the next section, a "cleaner" way to do the experiment 
via cavity locking finds no intermediate states exist between the first and second 
branches. 

In these plots, the point of onset of the second branch points reflect the atom 
numbers at which the second threshold appears for some detuning (not necessarily 
the detuning being considered). 
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Figure 5-25: Photon number n vs. effective atom number iV e ff for center of peak 
A = ±13.5MHz (cavity on resonant with atoms of most probable velocity). Circles 
(o): Data from positive-direction detuning scans with positive detunings, negative- 
direction scans with negative detunings. Crosses (x): Data from positive-direction 
detuning scans with negative detunings, negative-direction scans with positive detun- 
ings. Solid line: Rate equation model from Chap. EJ Dashed line: location of second 
threshold from corresponding quantum microlaser theory. 



130 




200 400 600 800 1000 



Figure 5-26: Photon number n vs. effective atom number N e g for A = ±12. 
Circles (o) and crosses (x) are as in Fig. 15-251 Line: rate equation model 
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Figure 5-27: Photon number n vs. effective atom number N e g for A = ±10. 
Circles (o) and crosses (x) are as in Fig. 15-251 
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Figure 5-28: Photon number n vs. effective atom number N e g for A = ±15. 
Circles (o) and crosses (x) are as in Fig. 15-281 
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Figure 5-29: Photon number n vs. effective atom number N e s for A = ±20.0MHz. 
Circles (o) and crosses (x) are as in Fig. 15-251 Line: rate equation model. Third 
threshold is apparent. 
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Figure 5-30: Photon number n vs. effective atom number N e g for A = ±7.5MHz. 
Circles (o) and crosses (x) are as in Fig. 15-251 Line: rate equation model 
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5.2.2 Disagreement with theory 

The data does not agree with theoretical curves in several ways. First, the theory 
underestimates the photon number in the case of cavity detuned significantly from 
resonance with the atoms. For example, see Fig. I5-3Q[ Fig. !5-27[ and Fig. 15-291 For a 
given density, the theory predicts a narrower first-branch resonance than is observed. 

Second, the peaks for the second and third branches occur at greater absolute 
detunings than for the initial branch. The third branch, in particular, exhibits a 
maximum at approximately ±28 MHz, more than twice the detuning of the first 
branch peak; the third threshold does not appear at all for absolute detunings less 
than about 17.5 MHz. From the model of Chapter 121 we would expect all branches to 
be approximately centered on zero detuning relative to atoms of the most probable 
velocity. 

In Sec. 16.51 we discuss possible reasons for the disagreement. 

5.2.3 Peak heights of detuning curves 

Suppose we consider the points of maximum photon number in the detuning curves 
as representing the condition of the cavity resonant with atoms of the most probable 
velocity. 

In Figure 13-311 we plot the peak in each detuning curve against the effective atom 
number. Transitions to the second and third branches occur very close to eachother. 
The agreement between theory and experiment is quite good except for the third 
branch, for which data points are considerably smaller than the theory values. Note 
also that the transitions to the third branch occur at a lower atom number the than 
predicted by the fully quantized theory. 

5.3 Cavity locking 

In order to study the system at a single detuning corresponding to resonance with the 
most probably velocity atoms, and to consider the effect of different initial conditions 
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Figure 5-31: Peak photon numbers of cavity scanning data vs. effective atom number. 
Solid line: rate equation solution for atom-cavity resonance (A = ±13.5MHz). Dotted 
line: location of thresholds by quantum theory. 
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of the cavity field more directly, we conducted experiments with the cavity locked on 
resonance. 

Cavity locking experiments used the following sequence: 

1. Lock cavity with a probe beam which has frequency corresponding to most 
probable velocity atoms, as described in Chapter |U 

2. Unlock cavity (remove probe beam and disengage cavity lock circuit) 

3. Move atom density modulator by a certain step size 

4. Measure atomic density by 553 nm probe and CCD (no 791 nm pump or probe) 

5. Lock cavity 

6. Unlock cavity 

7. Turn on the 791 nm pump beam, 

8. Trigger photon counter to measure microlaser emission for 100 msec ( "unseeded" 
data); download data to computer 

9. With the 791 nm pump remaining on, lock cavity again 

10. Unlock cavity 

11. Trigger photon counter to measure microlaser emission for 100 msec ("seeded" 
data); download data to computer 

12. Turn off 791 nm pump and repeat 

The repeated locking and unlocking of the cavity was necessary to refresh the 
locking in between other operations. A locking time of 300 msec was used. 

The frequency shift of the probe beam relative to the atomic line center was found 
to be +12.94 MHz. 

Note that for each step in density, two microlaser data samples are collected: (1) 
"unseeded" , in which the cavity is empty of photons when the pump beam is turned 
on, and (2) "seeded" in which the pump beam enters while the cavity has a very large 
(no ~ 10 12 ) number of photons due to the cavity locking beam. When the cavity lock 
is turned off, the microlaser field develops from this large-n initial condition. 
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Figure 5-32: Photon number n vs. effective atom number N e g for cavity locking 
experiment. Circles (o): "unseeded" data; Crosses (x): "seeded" data. Solid line: 
Rate equation model from Chap. El Dashed line: location of second threshold from 
microlaser quantum theory. 
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Figure 5-33: Detail of Fig. 15-321 in the initial threshold region 
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5.3.1 Cavity locking results 

The results of the cavity locking experiment are plotted in Figs. 15-321 and 15-331 The 
unseeded results agree with data from the cavity scanning experiment. The seeded 
results agree with the first branch for densities up to the onset of bistability; above 
this point the results agree with the second stable branch. 

Several outlying points in which the photon counts recorded were extremely high, 
n ~ 10 5 , were removed. These points were clearly artifacts because they occasionally 
appeared even when no atoms were present. They were most likely due to a problem 
in the sequencing electronics which on occasion allowed some of the locking probe 
beam to leak into the cavity during the microlaser experiment. 

5.4 Discussion 

The interpretation of the cavity locking results seems straightforward. The microlaser 
is becoming trapped in metastable states. The transition rate from a metastable to a 
stable solution is longer than the time of the experiment, which is about 100 msec and 
limited to that order of magnitude or smaller by the cavity re-locking requirement. 
Note that 100 msec 10 5 t cav is very long by the time scale of the cavity decay 
time. The estimate of a typical tunneling time between metastable and stable states 
is given by the Fokker-Planck analysis is fU\ : 

t g ~ t cav exp(aiV ex ) (5.1) 

where a is on the order of unity. In our experiment we have N ex w lOiVeff ~ 5000 in 
the region where the transition would be expected to occur spontaneously; this gives 
a transition time scale that is extraordinarily long. A calculation of Eq. 12.661 gives a 
transition rate W ~ 10~ 91 r c for N c r = 600 and zero detuning. By this calculation, 
spontaneous jumps will occur on time scales of 1 second or faster only for points 
extremely close (less than ~ 10~ 3 in relative distance) to points of marginal stability. 
Then as far as our experiments are concerned, the system would effectively behave 
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as if no spontaneous transitions at all occur. 

Some evidence suggests the situation may not be as simple as this. In the cavity 
scanning experiments, transitions from the first branch to and from the second often 
occurs with a abrupt bend in the curve; this occurs in both scan directions. This 
seems to suggest this transition occurs spontaneously rather than as a result of a 
branch becoming marginally stable and then disappearing. Note that the curve of 
solutions, as a function of detuning, must have a divergent slope at a marginally stable 
point, (see Fig. 15- 34 J) . It may be that other sources of noise (classical or quantum) 
drive fluctuations in the system faster than would be predicted by the Fokker-Planck 
analysis. We note that studies of bistability in the micromaser ^T] found transition 
rates much higher than predicted by a similar Fokker-Planck analysis. 

A similar analysis of hysteresis applies to the cavity scanning experiments. When 
plotted as photon number versus detuning for a fixed N e g, rate equation solutions 
appear as closed loops or curves as discussed in Chapter |21 

Transitions from the first branch to and from the second often occurs with a 
abrupt bend in the curve; this occurs in both scan directions. This seems to suggest 
this transition occurs spontaneously rather than as a result of a branch becoming 
marginally stable and then disappearing. Note that the curve of solutions, as a 
function of detuning or atom number, has a divergent slope at a marginally stable 
point. (Fig. EH) - 

5.5 Summary 

We have observed second and third thresholds in the microlaser. Theory and ex- 
periment are in remarkably good agreement for the resonant case. Highly history- 
dependent behavior is observed which suggests that the system becomes trapped in 
long-lived metastable states. 

Possible reasons for disagreement in the nonresonant case will be discussed further 
in the next chapter. 
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Figure 5-34: Arrows: transitions due to disappearance of a stable point. Dotted line: 
spontaneous transition from a metastable to stable state, rr-axis may represent atom 
number or detuning 



143 



144 



Chapter 6 



Discussion and Conclusions 



6.1 Multiple thresholds and the phase transition 
analogy 

The discontinuous nature of the higher thresholds and the existence of hysteresis sup- 
port the identification of these jumps as first-order transitions. Due to long metastable 
lifetimes, however, the observed jumps in the field do not occur at points predicted 
from the quantum microlaser theory. 

We may ask just how fruitful the phase transition analogy will be in this context. 
What insights from the study of many-particle systems can be applied to the micro- 
laser, or vice versa? Ideas from nonequilbrium statistical dynamics might be applied 
to the microlaser in describing the approach to steady-state and the transitions be- 
tween stable states. One theoretical paper [H7j suggests that under some variations a 
many-atom "mesolaser" with injected signal will exhibit not only multistability but 
temporal instabilities, period doubling, self-pulsing, quasiperiodicity, and chaos. 

We note that the microlaser is itself a many-particle system, but the particles are 
not strongly interacting; they may interact only through the common cavity field. 
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6.2 Many- atom effects 



The microlaser output in the resonant case agrees remarkably well with the rate 
equation solutions of Chapter |21 Both the quantum microlaser theory and the rate 
equation were derived by considering the influence of a single atoms on the field 
and assuming all atoms act independently. These experiments confirm the validity 
of this independent approximation. At first glance this may seem surprising, since 
the eigenvalues for the atom-cavity system with even two atoms are different from 
the single atom case. On the other hand, note that in general the photon number 
distribution in the cavity is strongly peaked (in fact, often more strongly peaked than 
a classical beam). The effect of the other atoms in the cavity can be neglected so long 
as their perturbation 5n is small compared to n. Moreover, the effect of the other 
atoms will tend to cancel eachother if the phases of their own Rabi oscillations are 
randomly or broadly distributed. 



6.3 Microlaser s, lasers, and randomness 

We now comment on the relation between the microlaser and a conventional laser. The 
microlaser's unique properties stem from its fully coherent atom-cavity interaction, 
which leads to an emission probability per atom (in the case of no broadening) of 

P e = sin 2 (Vn + lgtud (6.1) 

To model a conventional laser we simulate atomic decay by averaging over a transit 
time distribution weighted by (1/ti) exp(— t/r a ): 

1 r°° 1 / in + l)a 2 T 2 \ 

Pe = ~ dt int exp(-t/r«) sm 2 (V^TT<?tint) = - 1 , , ^ 2 2 ( 6 - 2 ) 
T a Jo 2 \1 + (n + l)g 2 r 2 J 

Solving the resulting gain-loss rate equation (Fig. I6-1J) gives a single solution for 
the photon number which increases linearly over threshold as a function of injection 
rate. Substituting this equation into the quantum microlaser theory gives steady-state 
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Figure 6-1: Gain-loss graphical solution for conventional laser as modeled in the text. 
Gain function follows saturation curve. 

photon statistics which are Poissonian far above threshold. 

From this perspective, the microlaser can be seen as a laser with an insufficient 
amount of noise or randomness to be a conventional laser. Conversely, a laser must 
possess a certain amount of noise to behave in the conventional way with monosta- 
bility and a single threshold. 

We note that any realistic microlaser will contain some averaging due to velocity 
distribution, coupling variation, etc. For a sufficiently large photon number, roughly 
n > / 5(gt- int )) 2 where 5(gt int ) is uncertainty in gt [nt the averaging will cause wash 
out oscillations in the gain function fFig J6-2|) to a value of one-half. Therefore the 
conventional laser is invariably recovered in the limit of large photon number. 

6.4 Value of quantum theory 

Our quantum theory predicts the steady-state photon number distribution in the mi- 
crolaser cavity. However, long metastable lifetimes prevent the true steady-state from 
being reached on an experimental time scale. Since the quantum theory otherwise 
agrees with the rate equation solutions, we might ask what value the quantum theory 
has in our case. One answer is that it is a more rigorous treatment than obtained by 
the rate equation, and therefore serves as a validation for it on theoretical grounds. 
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Figure 6-2: Gain-loss graphical solution for microlaser with some broadening. For 
sufficiently large n, gain line is flat and solution is similar to that of conventional 
laser. 

Another consideration is that the quantum theory also describes the statistics of the 
photon field. Experiments have been done to investigate the microlaser's photon 
statistics in the threshold region [HI] and an efforts to measure nonclassical statistics 
in the current microlaser are now under way. 

6.5 Detuning curves 

How do we explain the discrepancy between theory and experiment for the cavity 
nonresonant with the most probable velocity atoms? 

If we confine our attention to the resonance with a single traveling-wave compo- 
nent, the asymmetry in the detuning curves between positive and negative detunings 
(relative to the most probable velocities) is particularly puzzling. For a truly monov- 
elocity atom beam one would expect perfect symmetry between positive and negative 
detunings, by virtue of the symmetry in the Bloch equations when final state proba- 
bilities are measured. In our case we have a quite narrow beam: At> FW MM/^o ~ 17%, 
corresponding to a FWHM Doppler broadening of only 17% x 13.5MHz « 2.3MHz. 
And yet the center of the third branch (on the positive detuning side) is found at 
a detuning of about +25MHz relative to atomic line center, or +12MHz relative to 
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the resonance with the peak of the velocity distribution. It seems unlikely that the 
velocity distribution would lead to such a large effect. 

We now suggest a two ideas which are being explored in connection with the 
microlaser's detuning dependence. 



6.5.1 Influence of nonresonant traveling- wave field: atom in 
a bichromatic field 

The microlaser's traveling-wave (TW) interaction scheme requires a Doppler splitting 
larger than the widths of the microlaser resonances (i.e. no significant overlap between 
the two peaks). It was initially believed that if the two resonances are separated in 
this way, the interaction with the nonresonant TW component can be neglected. 
This may not be true: the nonresonant field may exert a significant influence on the 
resonant interaction even it would be directly responsible for an insignificant number 
of photons in the absence of the resonant interaction. 

To treat this problem exactly, we note that in the atom frame of reference, the 
cavity field consists of two components of frequency u C!W ± kv9. The equations of 
motion for an atom in a bichromatic field can be derived, analogous to Eqns. 12.151 
and EH 



2 2 



e -«02 e «(^O-^2)i 
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(6.3) 
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C a 



(6.4) 



where uj\ = u cav + kv9, u 2 = uj ciw + kv9 are the field frequencies, uj q is the atom 
resonance frequency, and 0i, 2 are the field phases. The solutions strongly depends 
on the choice of initial phase difference between the two fields. In a simulations, an 
averaging over phase difference must be performed. 
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6.5.2 Dispersive effects 



A related question is that of dispersive effects of the resonant and nonresonant atoms 
in the cavity. The susceptibility of a collection of iV atoms in a light field of frequency 
to (see jHH] ) is given by x — x' + *x"> with 



where N' is the atom density, Qr is the Rabi frequency of atoms in the field, 5l is 
the laser-atom detunig, and T a is the atom linewidth. 

This equation must be modified for the microlaser system, in which transit time 
is much shorter than atom lifetimes. A first approximation is to replace T a with t^ t . 

The refractive index of the atoms is 



n 



1 + x' » 1 + *72 (6-7) 



The shift in cavity frequency due to atom dispersion will be on the order of 
<5cav ~ x'/2c(j cav which may be as large as ~ 10 MHz for 1000 atoms in the cavity. The 
shift is a function of the Rabi frequency l^frig and atom-cavity detuning. In general 
there will be two such frequency shifts, one for each traveling- wave component, which 
add in effect. 

Frequency shifts from the dispersive effects may lead to another form of optical 
bistability, due to a nonlinear medium. The hysteresis observed in the microlaser 
may in fact be a combination of the microlaser's oscillatory gain function and "clas- 
sical" optical bistability effects. Unraveling the two forms of bistabilities will be an 
interesting challenge. 
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6.6 Future directions 



6.6.1 Theory of microlaser 

Clearly, the next step is to develop a more complete description of the microlaser, 
particularly for the nonresonant case. Additional data for high densities, and results 
from cavity locking experiments as a function of detuning as well as atom number, 
should provide important clues to the source of discrepancy between experiments and 
current theory. The possibility of many-atom effects playing a significant role cannot 
be ruled out, despite what the quantum trajectory results of Chapter El might suggest. 
In parallel, a number of other microlaser experiments are being considered. 

Second-order correlation measurements 

We have performed several experiments to measure the second-order correlation func- 
tion g^ 2 \t) of the microlaser emission 84 . Experiments have shown an reduction of 
amplitude of fluctuations, and a local maximum of the time scale of fluctuations, 
during the initial (second-order) threshold region. Our primary goal however has 
been to observe antibunching in the output field. Although the photon statistics in- 
side the cavity are highly squeezed in amplitude (with a Mandel Q parameter as low 
as -0.8), the resulting nonclassical effect in the output light is limited by the rela- 
tively large number of photons in the cavity. The second-order correlation function 
at zero time delay is related to the average cavity photon number and Q parameter 
by gW(0) = 1 + Q/(n); for (n) « 500 we expect gW(0) - 1 « -0.002. Classical 
fluctuations and other systematic effects may obscure this small effect. 

Internal statistics: Cavity dumping experiment 

One method of more directly investigating photon statisics of the cavity field is via a 
"cavity dumping" scheme. The cavity-atom interaction could suddenly be disabled, 
by detuning or by optical mean, after which the internal photons are collected as they 
exit the cavity. By repeating this process many times, a photon count distribution 
P^f CD ) is obtained which is related to the cavity photon number distribution p^( PND ) 
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by the partition formula 



p(PCD) _ 



m / j 

n=m 



oo 




f(l-r,) n ?( PND ), 



(6.8) 



where 77 is the overall counting efficiency. 
Microlaser lineshape measurement 

Lineshape measurements in the many atom regime is a potentially rich field of study. 
Number-phase uncertainty requires sub-Poissonian statistics to be associated with 
increased phase diffusion, and thus a larger linewidth. For this reason we expect to 
observe dramatic changes in microlaser linewidth as the number of atoms and other 
parameters are varied. Such behavior has been predicted in 03] in the context of the 
micromaser, but has yet to be observed, either directly or indirectly. 

6.7 Summary 

We have developed a laser system in which atom-field interaction is truly coherent. 
Rabi oscillations of the atoms lead to multistability in the laser field intensity, and 
rapid jumps occur between different stable points. Theory and experiment are in 
good agreement for the case of the cavity resonant with atoms of the most probable 
velocity, but long lifetime of metastable states prevent verification of true steady-state 
transition points. Hysteresis is observed as a function of atom-cavity detuning and 
initial cavity photon number. The detuning behavior of the microlaser is not yet 
well-understood; interactions with the nonresonant field and dispersive effects of the 
resonant and nonresonant atoms are being investigated. 

Quantum trajectory simulations support the application of single-atom theory to 
the many-atom case. In addition, simulations predict an increase in photon number 
variance proportional to the product of the interaction time and cavity decay rate, 
for parameters such that the single-atom theory predicts sub-Poisson statistics. 
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Appendix A 



Microlaser theory calculations 



These programs were written and executed using Mathematica 3.0 and 4.0 by Wolfram 
Software. 

A.l Micromaser theory 

(* All distances in microns *) 
(* All times in microseconds *) 
(* All frequencies in 10~6 /s *) 

us = 1; 

mhz = 2. Pi; 

lambda = 0.7911; 
wm = 41.1 ; 

vO = 815; (* most probable velocity, m/s *) 
deltapeak = +13.5 mhz; 

tiltangle = (deltapeak / (2Pi)) * lambda / vO; 
gammacav = 1/0.88; 
v = 815. ; 

tint = Sqrt[Pi] * wm / v; 
deltaVOverV = 0.165; 

gcouplingO = 2*Pi*192.9 * 10~-3; 
gcoupling = gcouplingO ; 

deltacavlowerlimit = 0; 
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deltacavupperlimit = 50 mhz; 
sqnupperlimit = 150; 



doMicromaserCalc := 

(pOverC = Table [0, {k, 0, kmax}] ; 
p = Table [0, {k, 0, kmax}]; 
p0verC[[0 + 1]] = 1; 
sumpOverC = 1 ; 
Do[ 

pOverC[[k + 1]] = pOverC[[k]] * Nex * 
betaapprox [Sqrt [k] , delta] /k ; 

sumpOverC = sumpOverC + pOverC [ [k + 1] ] , 
{k, 1 , kmax}] ; 

Do[p[[k]] = pOverC[[k]] / sumpOverC, {k, 1, kmax + 1}]; 
nAvg = Sum[k * p[[k + 1]], {k, 0, kmax}]; 

n2Avg = Sum[k~2 * p[[k + 1]], {k, 0, kmax}]; 
mandelQ = (n2Avg - nAvg" 2 - nAvg) /nAvg; 

) ; 

getnQ := (Nex = natoms / (0.10595 * 1); 

kmax = Floor [Nex + 2]; doMicromaserCalc; Return [{nAvg , mandelQ}] ) ; 



A. 2 Calculation of gain function 

This code includes Dichromatic effects (set fi 2 = for ordinary monochromatic theory) 



Deltal =10; (* relative to atoms *) 

Delta2 = 37 mhz; (* relative to atoms *) 
tstart = -2* tint; 
tstop = +2 * tint; 

dosoln := NDSolve [{ca' [ 

t] == ((I *omegaRl [t] /2) E~(-I *phil) E~(I* Deltal* t) + (I * 
omegaR2[t]/2) E~(-I *phi2) E~(I Delta2 * t))* cb[t], 
cb'[t] == ((I *omegaRl [t] /2) E~(I *phil) E~(- I* Deltal* t) + (I * 
omegaR2[t]/2) E~(I *phi2) E~(- I Delta2 * t))*ca[t], 
ca [tstart] == 1., cb [tstart] == 0.}, {ca, cb}, {t, tstart, tstop}, 
MaxSteps -> 10000, Start ingStepSize -> Automatic, 
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MaxStepSize -> tint/10, WorkingPrecision -> 10] ; 
(* Average over relative phases! *) 

numphases = 10 . ; 

phase [i_] := 2. Pi * (i - 1) /numphases ; 
solarray = Table [i, {i, 1, numphases}]; 

dophaseavg := 
( 

omegaO = 2* Sqrt [n + 1] * gcoupling; 
omegaRl [t_] := omegaO * Exp[-(v t/wm)~2] ; 
omegaR2[t_] := omegaO * Exp[-(v t/wm)~2] ; 

Do[ 
( 

phil = 0*phase [i] ; 
phi2 = phase [i] ; 

solarray [ [i] ] = dosoln;), {i, 1, numphases}]; 

caarray = Table [(ca /. solarray [ [i] ])[ [1] ] , {i, 1, numphases}]; 

(* cbarray = Table [(cb /. solarray [ [i] ])[ [1] ] , {i, 1, numphases}]; *) 

ca2avg2[t_] := Sum [(Abs [caarray [ [i] ] [t] ] "2) /numphases, {i, 1, numphases}]; 
(* cb2avg2[t_] := 

Sum [(Abs [cbarray [ [i] ] [t] ] "2) /numphases , {i, 1, numphases}]; *) 
(* ca2avg[t] = 

Functionlnterpolation [ 

Sum [(Abs [caarray [ [i] ] [t] ] "2) /numphases , {i, 1, numphases}], {t, 
tstart, tstop}] ; *) (* 

cb2avg = 

Functionlnterpolation [ 

Sum [(Abs [cbarray [ [i] ] [t] ] "2) /numphases , {i, 1, numphases}], {t, 
tstart, tstop}]; *) 
Return [1 - ca2avg2 [tstop] ] ; 

) ; 

gaussian[x_, xa_, sigma_] := 

l/(Sqrt[2. Pi]*sigma) E~(-(x - xa) "2/(2 . *sigma~2) ) ; 
veldist [vel_] := gaussian[vel, vO, v0*deltaV0verV / 2.355]; 

numvels =11; (* odd number is better *) 
vlowerlimit = vO (1 - 1 . 2*deltaV0verV) ; 
vupperlimit = vO (1 + 1 . 2*deltaV0verV) ; 
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vsample [i_] : = 

N [vlowerlimit + (i - 1) * (vupperlimit - vlowerlimit) / (numvels - 1)]; 

deltaof v [vel_] := (vel/vO - 1 . ) *deltapeak; 
vweight[i_] : = 

veldist [vsample [i] ] * (vupperlimit - vlowerlimit) / (numvels - 1); 
vlist = Table [{vsample [i] , vweight [i]}, {i, 1, numvels}]; 

numgs = 4; 
gvariation = 0.09; 
gupperlimit = . 99*gcoupling0; 

glowerlimit = 1.01 * gcouplingO* (1 - gvariation); 
gdist[g_] : = 

UnitStep[l - g/gcouplingO] *UnitStep [g/gcouplingO - (1 - gvariation)]* 
Sqrt [Log [gcouplingO/g] ] ; 

gsample[i_] := 

N [glowerlimit + (i - 1) * (gupperlimit - glowerlimit) / (numgs - 1)]; 
gweight2[i_] := gdist [gsample [i] ] * (gupperlimit - glowerlimit) /(numgs - 1) ; 
gnorm = Sum[gweight2 [i] , {i, 1, numgs}]; 
gweight[i_] := gweight2 [i] /gnorm; 

glist = Table [{gsample [i] , gweight[i]}, {i, 1, numgs}]; 

(* betal : averaged over phase, velocity = detuning *) 

dobetal := 

Sum [ vweight [i] * 

( v = vsample [i] ; 

Deltal = deltacav - (v/vO) * deltapeak; 

Delta2 = -(deltacav + (v/vO) *deltapeak ); 
dophaseavg 
) , {i , 1 , numvels}] ; 

(* beta2 : averaged over phase, velocity = det, gcoupling *) 

dobeta2 := 

Sum[gweight [j] * 

( gcoupling = gsample [j]; 
Sum [ vweight [i] * 

( v = vsample [i] ; 
Deltal = deltacav - (v/vO) * deltapeak; 

Delta2 = -(deltacav + (v/vO) *deltapeak ); 
dophaseavg 
), {i, 1, numvels}]) 
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, { j , 1 , numgs}] ; 
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Appendix B 



Quantum trajectory program 



B.l Main program: qts.c 



#include <stdio.h> 
#include <math.h> 

#include <stdlib.h> /* for exit(), rand() , and mallocO */ 
#include <time.h>/* for clock() and time() */ 

/* Version 9.5 (CFY) : include comments in input files */ 
/* Fixed bugs in counting number of atoms, average number of photons */ 

/* This is version 9, single precision version. This version 
includes atom-cavity detuning effect. For this _C[] becomes 
_C[] (real) and _D[] (imaginary). */ 

/* adaptive step size control: dt is not constant */ 

/* Unused defs */ 

/* #define USE_MSDOS_NAMES */ 

/* #define REUSE_C_ARRAY 

/* initialize _C[] with the last saved _C[] array */ 
/* #define CLICK_INTERRUPT 

/* sense mouse click */ 
/*#define X_GAUSSIAN_BASED_ON_RABI_OSC 1 */ 

/* 1 means that <N> calculation is based on Rabi oscillation */ 
/* argument, i.e., averaging of g(x) , whereas means that it is */ 
/* based on averaging of g~2. */ 
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#define SUMMARY_F I LENAME_DEFAULT "qts.data" 

#define USE_RK4_ALG0RITHM /* Use Runge-Kutta 4th order */ 

#define DRAND48.EXIST 

/* reproducible random numbers. */ 
/* displaying dots */ 



#define RANDOM.INIT 
#define SHOW.PROGRESS 



#define LINUX 

#define VERSION 9 

#define REVISION 1 

#define A_INFO_ITEMS 7 

#define MAX_NUM_ATOM 20 

#define MAX_NUM_INPUT_FILES 64 

#define MAX.MEMORY two_to(25) 



/* 22-> 4 MB, 23-> 8MB, 24->16MB */ 



#define D_STATE_DECAY_CORRECTION 1.34 

/* D-state decay ratein 1P1-1S0 decay rate 



amj 1.0 -> 1.37 */ 



#define WAIT_PORTION_BEFORE_AVERAGING 0.0 

/* amj changed it to .0 instead of . 1 : */ 
/* . 1 = first 10% not used for averaging */ 



/* same thing as above, but active in this version (CFY) */ 

#define VERY.SMALL l.e-30 
#define SQRT.CORRECTED 

/* Platform specific redefinition */ 

#ifdef DRAND48.EXIST 
#define drand() drand48() 
#define srand(x) srand48(x) 
#else 

#define drand() (rand() / (double) RAND.MAX) 

/* for Mac and PC */ 

#endif 

#ifdef SQRT.CORRECTED 
#define sqrt(x) sqrt(fabs(x)) 
#endif 

#define sqr(x) ((x)*(x)) 

#ifdef LINUX 
#define pi PI 
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#else 

#define PI pi 
#endif 

#define PI 3.1415926536 

/* output/input file names */ 

#define QT.LOG "test. log" 
#ifdef USE_MSDOS_NAMES 
#define INDEX_FILE "test.idx" 
#define INPUT.FILE "test.inp" 
#define STATUS.FILE "test.sta" 
#define OUTPUT.FILE "test.out" 
#define PHOTON.FILE "test.pho" 

#else 

#define INDEX_FILE "test. index" 
#define INPUT_FILE "test. input" 

#define STATUS.FILE "test . status" 

#define OUTPUT.FILE "test.output" 
#define PHOTON.FILE "test.photon" 
#endif 

/* global variables accessed from other routines */ 
double *_C, *_D, *_C_1, *_D_1; 

#ifdef USE_RK4_ALG0RITHM 

double *_C_2, *_C_3, *_D_2, *_D_3; 

#endif 

/* C=real, D=imaginary part */ 

double _u_th, _w0, _Gc_over_2pi , _g_over_2pi, _Ga_over_2g; 

int _max_atoms, _max_photon_size; 

double _A_inf o [MAX_NUM_ATOM] [A_INFO_ITEMS] ; 

int _N_atom; /* current number of atoms in the cavity */ 

double _t_int, _Delta_t, _beam_radius , _delta_v, _v0, _detuning=0 . ; 

/* amj add; _v0 is peak of selected velocity in most 

/* probable velocity unit */ 

double *SIN, *EXP; 

double *_Pn_dist, *_Pn_dist_mean, *_Pn_dist_mean_beginning; 

/* Note: Pn_dist_mean is actually P(n) * (averaging time / dt) */ 

/* until end of sim. */ 
double .population [MAX_NUM_ ATOM] ; 



/* for INPUT.ONLY option */ 
/* contains _C[] array, etc */ 
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double _mean_shif t_over_2g, _CA_detuning; 
double **_bits, *_detuning_sum, _omega_Rabi; 
double g_multiplier ; 
long _max_k; 

/* #ifdef SPONT_EMISS_IN_PUMP 
double _pump_waist=l .0; 
#endif */ 

/* User-specified parameters */ 

int INCLUDE_ATOM_DECAY ; /* atomic spontaneous emission */ 
int INCLUDE_D_STATE_DECAY; /* D-state decay loss included */ 
int USE_LAST_C_ARRAY=0 ; 

int STANDING. WAVE; 

/* include sin kz dependence in g */ 

int GAUSSIAN_X=1; /* amj jul 27, 98 */ 

/* assume a gaussian exp(-x~2/w0~2) dependence */ 
/* along x axis. */ 

double INTERACTI0N_LENGTH=2. ; 

/* atom-field interaction length. */ 

/* should be larger than 2 in _w0 unit */ 

int GAUSSIAN_Y=1; 

/* assume a gaussian exp(-y~2/w0~2) dependence */ 
/* along y axis */ 

int CIRCULAR.BEAM; 

/* assume a circular atomic beam with diameter of D */ 
double REC_HALF_HEIGHT ; 

/* half height of the beam in _w0, width is equal to D */ 
int M0N0_VEL0CITY=0; 

/* mono-energetic atomic beam with v=u_th */ 

int USE_MAXWELL_BOLTZMANN ; 

int DEBUG_LEVEL=0 ; 

/* 0=only summary l=immediate output 2=debug info 
/* 3=extra info */ 
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int PUMPING.SIMULATED; 

/* pump field effect during evolution included */ 



double 



PUMPING.LENGTH; 
/* should be large enough to include 2*pump waist */ 



double 



PUMP.WAIST; 

/* waist of pump field assuming gaussian */ 



double 



PUMP.POS ; 

/* pump center to cavity center distance. */ 
/* should be positive. */ 

/* integration from -PUMPING.LENGHT/2 to */ 

/* PUMP_P0S-INTERACTI0N_LENGTH/2 with the pump */ 

/* center is located at 0. */ 



int PUMPING_EFFECT_ON=0; 



int 



n_seed=0, 

f irst_atom_f lag=0 ; 



char *inputline, *inputf ilename, *indexf ilename , *statusf ilename , 



*outputf ilename , *photonf ilename , *pnf ilename, *summaryf ilename ; 



/* Function declarations */ 

int getline(FILE *ff, char *s, int lim) ; 
void testroutine(void) ; 
double get_arrival_time(void) ; 

/* double get_velocity (void) ; */ 

/* double get_arrival_time_old(void) ; */ 

main (int argc, char **argv) 
{ 

double Ga_over_2pi, lambda, rO, finesse, L; 

double N_eff[100], u_th[100] , zO, fsr, f, t, dt, Ti, Ti2, Ti2new; 
double Tl, T2, time_elapsed; 



double 

int 

int 

double 



N_ex, theta_pump; 
VARIABLE.FINESSE = 0; 
ADJUST_N_EFF = 0; 
new_N_ef f ; 



double 
int 



WAIT.FRACTION =0.3; 
WAIT_FRACTION_TENTHS ; 
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long ii, C_size, total_array_size, max_k, P_atom [MAX_NUM_ATOM] ; 

int i, j, num_N_eff, t ot_ at om_ count , m, atom_count, 

atom_count_since_Ti2, photon_count , too_many_atom_count , 
num_u_th, mm, injected_state, num_excited_exiting; 

int atom_decay_to_d_count=0; /* amj */ 

time_t tp; /* for initializing srandO */ 

double D, next_entry_time , t_int, delta_t, n_mean, n2_mean, 

n_aver, n2_aver, n_sum, n2_sum, X, n_meanl, n_mean2; 

/* amj add n2_mean and change dn_sum to n2_sum and */ 

/* dn_aver to n2_aver */ 
double P_sum, P_N_mean, P_N_mode, emission_prob; 

FILE *fp; 

int MAKE_INPUT_ONLY=0 , file.id; 

char output_name [80] , output_time [80] , photon_f ile_name [80] ; 

char *inputfile[MAX_NUM_INPUT_FILES] ; 

int NUM_INPUT_FILES=0, K, GENERATE_PH0T0N_FILE=0 ; 

double num_excited_injected_atoms , excited_state_prob; 

/* expectation value */ 

int spont_count=0, d_state_decay_count=0, stop_now=0; 

double tilt_angle, detuning [100] , mean_shift; 

int num_detuning , mmm, nonzero_tilt ; 

double vacuum_time, ef f ective_V_to_apparent_V, 

big.dt, dtO, dtO.init, Mandel.Q; 
long num_dt_steps ; 

double num_dt_aver_steps=0 . ; /*amj number of dt steps */ 
/* during which average of n_aver and n2_aver is taken */ 
double num_dt_aver_steps_last=0 . ; /* placeholder */ 
double n_final[100] , Q_f inal [100] , N_f inal [100] , Gc_f inal [100] ; 

double arrt, p_poisson; 

/* function declarations */ 

void evolve_wave_f unction (double , double); 

long two_to(int); 

void def ine_sin_exp() ; 

void read_status_f ile() ; 

void write_status_f ile () ; 
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double inj ect_new_atom (double) ; /* returns excited state prob. 
int cavity_decay_occurs (double, double* , double*) ; 

/* amj add , double* */ 

int how_many_atoms_decay(double) ; 
int how_many_d_state_decay (double) ; 

double final_n_avg; 
double f inal_n2_avg; 
double f inal_normalize; 

#ifdef CLICK. INTERRUPT 
int key_pressed() ; 
#endif 

/* #ifndef LINUX */ /* amj */ 
float volume_ratio () ; 
/* argc=ccommand(&argv) ; */ /* amj */ 
/* #endif */ /* amj */ 

/* make sine table */ 
/*def ine_sin_exp() ;*/ 

/* Memory allocation for filename strings */ 

inputline = malloc(100 * sizeof (char) ) ; 
inputf ilename = malloc(15 * sizeof (char) ) ; 
/* inputline = malloc(20 * sizeof (char) ) ; */ 
inputf ilename = malloc(20 * sizeof (char) ) ; 
indexf ilename = malloc(20 * sizeof (char) ) ; 
statusf ilename = malloc(20 * sizeof (char) ) ; 
output filename = malloc(20 * sizeof (char) ) ; 
phot onf ilename = malloc(20 * sizeof (char) ) ; 
summaryf ilename = malloc(20 * sizeof (char) ) ; 
pnf ilename = malloc(20 * sizeof (char) ) ; 



/* printf("(0) Top hat or (1) Gaussian in x direction?\n") ; 
scanf 07.d", &GAUSSIAN.X) ; */ 

/* Get an input file list */ 

for (i=0; i<MAX_NUM_INPUT_FILES; ++i) 
inputf ile [i] =malloc (80) ; 
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fprintf (stderr , "Enter input file name(s) \n"); 

f printf (stderr , "(max number of names < %d and '&' should be the 
last item in the list.)\n\ MAX_NUM_INPUT_FILES) ; 
fprintf (stderr, without the list "); 

fprintf (stderr , "will make a new input file.\n"); 
fprintf (stderr , "Enter $ for test routine . \n") ; 

i=-l; 

NUM_INPUT_FILES=0; 
do { 

i++; 

scanf ("°/„s" , inputf ile [i] ) ; 

if (*inputf ile [i] == '$') testroutine () ; 
} while ( *inputfile[i] != '&'); 

NUM_INPUT_FILES=i; 

if (i==0) MAKE_INPUT_0NLY=1 ; 

if (MAKE_ I NPUT_ ONLY ) 

printf ( " input only\n" ) ; 

else 

{ for (i=0; i<NUM_INPUT_FILES ; ++i ) 

{ fprintf (stderr , "batch %d: %s\n" , i, inputf ile [i] ) ; 
f p=f open(inputf ile [i] , "r"); 
if (f p==NULL) 
{ printf ("file not found, exiting\n"); 
/* fclose(fp); */ 
exit (1) ; 

} 

f close(fp) ; 

} 

} 



if (MAKE_INPUT_ONLY) 
{ 

printf ("Enter new input filename (max 15 chars) :\n"); 
scanf ("°/„s", inputf ilename) ; 

f p=f open(inputf ilename , "w"); 

/* user input */ 

fprintf (stderr , "Microlaser quantum trajectory simulation"); 
fprintf (stderr, "version °/ d\n" , VERSION); 
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fprintf (stderr , "Enter Ga/2pi in kHz:\n"); 
scanf ("°/,lf " , &Ga_over_2pi) ; 
fprintf (fp, "°/ lf # Ga_over_2pi\n" , Ga_over_2pi) ; 

fprintf (stderr , "Enter wavelength in nanometer : \n" ) ; 

scanf C7.1f\ felambda); 

fprintf (fp, "°/„lf # lambda\n", lambda); 

fprintf (stderr , "Enter radius of curvature in cm:\n"); 

scanf ('7„lf\ &r0); 

fprintf (fp, "°/„lf # rO\n", rO) ; 

fprintf (stderr , 

"Enter finesse in million, negative for variable : \n") ; 
scanf ("%lf", fefinesse); 
fprintf (fp, "°/„lf # finesse\n", finesse); 

fprintf (stderr , "Enter cavity length in cm:\n"); 

scanf C7.1f\ &L); 

fprintf (fp, "°/„lf # L\n", L); 

fprintf (stderr , "Enter max number of photons +1 and "); 
fprintf (stderr , "atoms in the cavity:\n"); 

scanf ("°/ d °/ d" , &_max_photon_size, &_max_atoms) ; 
fprintf (fp, "°/„d °/„d # _max_photon_size _max_atoms\n" , 
_max_photon_size , _max_atoms) ; 

fprintf (stderr , "Enter (initial) size of integ. step in l/2g:\n"); 

scanf ("%lf", &dtO); 
fprintf (fp, "°/„lf # dtO\n", dtO) ; 

fprintf (stderr , "Enter total # atoms to inject into cavity:\n"); 

scanf ("°/ d" , &t ot _ at om_ count) ; 
fprintf (fp, "°/„d # tot_atom_count\n" , tot_atom_count) ; 

/* control parameters */ 

fprintf (stderr, "Enter a value for DEBUG_LEVEL(0 . .3) :\n") ; 
scanf ( "%d" , &DEBUG.LEVEL) ; 

fprintf (fp, "°/„d # DEBUG_LEVEL\n", DEBUG.LEVEL) ; 

#ifdef REUSE_C_ARRAY /* Reusing _C[] values */ 
fprintf (stderr, 

"Initialize _C[] to those from previous loop/run? \n"); 
fprintf (stderr, " (yes=l/no=0) :\n") ; 
scanf ("°/ d" ,&USE_LAST_C_ARRAY) ; 
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fprintf (fp, "°/„d # USE_LAST_C_ARRAY\n" ,USE_LAST_C_ARRAY) ; 
#endif 

/* damping processes */ 

fprintf (stderr , "Include atomic spontaneous emission? "); 
f printf (stderr , " (yes=l/no=0) : \n" ) ; 
scanf ( l,0 /„d" , &INCLUDE_ATOM_DECAY) ; 

fprintf (fp, "°/„d # INCLUDE_ATOM_DECAY\n" , INCLUDE_ATOM_DECAY) ; 

fprintf (stderr, "Include D-state branching? (yes=l/no=0) : \n") ; 
scanf ( "%d" , &INCLUDE_D_STATE_DECAY) ; 

fprintf (fp, "°/„d # INCLUDE_D_STATE_DECAY\n" , INCLUDE_D_STATE_DECAY) ; 

/* cavity mode parameters */ 
fprintf (stderr, 

"Choose (0) traveling wave, (1) standing wave mode:\n"); 
scanf C7„d" , &STANDING.WAVE) ; 
fprintf (fp, "%d # STANDING_WAVE\n" , STANDING.WAVE) ; 

if ( ! STANDING.WAVE) 

{ /* traveling-wave must have a nonzero tilt */ 
fprintf (stderr, 

"Enter atomic beam tilt angle(theta) in mrad:\n"); 
scanf ("°/ lf", &tilt_angle) ; 

fprintf (fp, "7„lf # tilt_angle\n" , tilt.angle) ; 
} 

else 

tilt_angle=0 . ; /* no tilt for standing-wave mode */ 
fprintf (stderr, 

"Enter number of atom-cavity detuning values (max of 100) :\n"); 
scanf ( "%d" , &num_detuning) ; 

fprintf (fp, "°/„d # num_detuning\n" , num_detuning) ; 

if ( ( ! STANDING.WAVE) && (f abs (tilt .angle) >1 . e-10) ) 
{ 

/* traveling-wave & nonzero tilt angle */ 

nonzero_tilt=l ; 
fprintf (stderr, 

"Enter detuning values in (u/lambda*theta) unit:\n"); 

} 

else 

{ 

nonzero_tilt=0 ; 
fprintf (stderr , "Enter detuning values in MHz:\n"); 
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} 



for (i=0; i<num_detuning; ++i) { 
scanf ("%lf ", fedetuning [i] ) ; 

fprintf (fp, "°/„lf # detuning [°/„d] \n" , detuning [i] , i) ; 
} 

fprintf (stderr , 

"Enter atom-cavity interaction \"full\" length in wO:\n"); 
scanf ("%lf " , &INTERACTION.LENGTH) ; 

fprintf (fp, "'/.If # INTERACTION_LENGTH\n" , INTERACTION.LENGTH) ; 

fprintf (stderr , "Choose (0) rectangular beam, (1) circular beam:\n"); 
scanf ("%d" , &CIRCULAR.BEAM) ; 

fprintf (fp, "°/od # CIRCULAR_BEAM\n" , CIRCULAR_BEAM) ; 

if (CIRCULAR.BEAM) 

{ 

fprintf (stderr, 

"Enter diameter of atomic beam D in um:\n"); 
scanf ("%lf", &D); 

fprintf (fp, "°/„lf # D diameter of beam\n" , D) ; 

} 

else 

{ 

fprintf (stderr, 

"Enter half height of atomic beam in w0:\n"); 
scanf ("°/„lf" , &REC_HALF_HEIGHT) ; 

fprintf (fp, "°/ lf # REC_HALF_HEIGHT\n" , REC_HALF_HEIGHT) ; } 

fprintf (stderr , "Ratio of effective volume to apparent volume?\n"); 
fprintf (stderr , "On Mac it is automatically calculated : \n" ) ; 
scanf ("°/ lf " , &ef f ective_V_to_apparent_V) ; 
fprintf (fp, "°/„lf # ef f ective_V_to_apparent_V\n" , 
ef f ective_V_to_apparent_V) ; 

/* End of cavity input */ 



/* N_eff[] input */ 

fprintf (stderr , "Enter number of <N> values (max of 100) :\n"); 
scanf ("°/od", &num_N_ef f ) ; 

fprintf (fp, "°/ d # num_N_ef f \n" , num_N_ef f ) ; 

fprintf (stderr , "Enter <N> values :\n"); 
for (i=0; i<num_N_eff; ++i) 
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{ scanf ("°/„lf ", &N_eff[i]); 
fprintf(fp, "°/ lf # N_ef f [°/„d] \n" , N_ef f [i] , i) ; 
> 

f printf (stderr , 

"Enter number of most probable velocity values (max of 100) :\n"); 
scanf ("°/„d" , &num_u_th) ; 

fprintf(fp, "°/ d # num_u_th\n", num_u_th) ; 

f printf (stderr , "Enter most probable velocities (u) in m/s:\n"); 

for (i=0; i<num_u_th; ++i) 

{ 

scanf ('"/.If", &u_th[i]); 

fprintf(fp, "°/olf # u_th[°/„d]\n\ u_th[i] , i) ; 

} 

f printf (stderr, 

"Choose (0) velocity distribution, (1) mono velocity : \n") ; 
scanf ("°/ d" , &M0N0.VEL0CITY) ; 

fprintf(fp, '7„d # M0N0_VEL0CITY\n" , M0N0_VEL0CITY) ; 

/* For non-Maxwell Boltzmann velocity distribution */ 
if ( ! (MONO.VELOCITY) ) 
{ 

f printf (stderr , "Which distribution\n") ; 

f printf (stderr , "(0) M-B convoluted with Lorentzian\n") ; 

f printf (stderr, "(1) Maxwell-Boltzmann? : \n") ; 

scanf ("°/ d" , &USE_MAXWELL_BOLTZMANN) ; 
fprintf(fp, "°/ d # USE_MAXWELL_BOLTZMANN\n" , 
USE_MAXWELL_BOLTZMANN) ; 

if ( ! (USE_MAXWELL_BOLTZMANN) ) 

{ f printf (stderr , "Enter full width of velocity "); 

f printf (stderr , "spread in most probable velocity unit\n"); 
/* scanf ("%lf °/ lf", &_delta_v, &_v0) ; */ 
scanf C7.1f", &_delta_v); 

f printf (stderr , "Enter peak of the selected velocity"); 
f printf (stderr, " in most probable velocity unit:\n"); 
scanf C7.1f\ &_v0); 

fprintf(fp, "°/„lf # _delta_v\n", _delta_v) ; 
fprintf(fp, "°/„lf # _v0\n", _v0) ; 
} 

} 
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fprintf (stderr , "Generate a record for photons emitted? 

(yes=l, 0=no)\n"); 

scanf ("%d" , &GENERATE.PHOTON.FILE) ; 

fprintf (fp, "°/ d # GENERATE_PHOTON_FILE\n" , GENERATE_PHOTON_FILE) ; 
/* 

fprintf (stderr , "Enter a & b where <n>a*<N>~b:\n") ; 
scanf ("%lf %lf", fea.coeff, feb.coeff ) ; 
fprintf (fp, "°/„lf °/ lf\n", a.coeff, b.coef f ) ; 
*/ 

/* pumping simulation, if required */ 

fprintf (stderr , "Include pump field effect on atom-cavity "); 
fprintf (stderr , "interaction? (yes=l/no=0) : \n") ; 
scanf ("%d" , &PUMPING.SIMULATED) ; 

fprintf (fp, "%d # PUMPING_SIMULATED\n", PUMPING.SIMULATED) ; 

if (PUMPING.SIMULATED) 

{ fprintf (stderr , "Pump field effect included: \n") ; 

fprintf (stderr, 
"Enter pumping length, pump waist in wO unit:\n"); 

scanf ("°/ lf %lf", &PUMPING.LENGTH, &PUMP_WAIST) ; 

fprintf (stderr, 

"Enter position of the pump in wO unit (+ if prior to the cavity) : \n") ; 
/* PUMP_P0S is positive if the pump is place before the cavity */ 
scanf ("°/ lf" , &PUMP.P0S) ; 
fprintf (fp, 

"°/olf %lf '/.If # PUMPING.LENGTH, PUMP.WAIST, PUMP_P0S\n" , 
PUMPING.LENGTH, PUMP.WAIST, PUMP.POS) ; 
> 

fprintf (stderr , "Enter filename for summary data, blank for ,0 / s':\n", 
SUMMARY_F I LENAME_DEFAULT ) ; 
scanf ("°/ s" , summaryf ilename) ; 

fprintf (fp, "°/ s # summary data f ilename \n", summaryf ilename) ; 
f close(fp) ; 

fprintf (stderr, "*** End of Input ***\n") ; 
exit(l) ; 
> 

/* From now on it is "not" MAKE. INPUT. ONLY */ 
/* Grand Loop for input files */ 
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for (K=0; K<NUM_INPUT_FILES; ++K) 

{ /****************** GRAND LOOP *******************/ 



/* Read from inputfile[K] */ 

if ((fp=fopen(inputfile[K] , "r"))==NULL) 

{ fprintf (stderr , "*** %s cannot be found. \n" , inputf ile [K] ) ; 
exit(l) ; 
} 

getline(fp, inputline, 100); 

sscanf (input line , "°/ lf", &Ga_over_2pi) ; 

/* fscanf(fp, " /„lf", &Ga_over_2pi) ; */ 

/* printf ("Ga_over_2pi = °/„f\n" , Ga_over_2pi) ; */ 

getline(fp, inputline, 100); 
sscanf (inputline , "°/,lf", felambda) ; 

getline(fp, inputline, 100); 
sscanf (inputline, "°/ lf", &r0) ; 

getline(fp, inputline, 100); 
sscanf (inputline , "°/ lf", fefinesse); 

if (f inesse<0) 

{ finesse = -finesse; 
VARIABLE.FINESSE = 1; 

} 

getline(fp, inputline, 100); 
sscanf (inputline , "°/ lf", &L) ; 

getline(fp, inputline, 100); 

sscanf (inputline , "°/ d °/ d" , &_max_photon_size , &_max_atoms) ; 

getline(fp, inputline, 100); 

sscanf (inputline , "°/ lf", &dt0_init) ; 

getline(fp, inputline, 100); 

sscanf (inputline , "°/„d" , &tot_atom_count) ; 

/* printf ("tot_atom_count = °/„d\n" , tot_atom_count) ; */ 

getline(fp, inputline, 100); 

sscanf (inputline, ,|0 /,d" , &DEBUG.LEVEL) ; 
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#ifdef REUSE_C_ARRAY 

getline(fp, inputline, 100); 

sscanf (input line, ,,0 /.d" ,&USE_LAST_C_ARRAY) ; 
#endif 

getline(fp, inputline, 100); 

sscanf (inputline, "°/ d" , &INCLUDE_ATOM_DECAY) ; 

getline(fp, inputline, 100); 

sscanf (inputline, "Id", &INCLUDE_D_STATE_DECAY) ; 

getline(fp, inputline, 100); 

sscanf (inputline, "%d" , &STANDING.WAVE) ; 

/* printf ("STANDING_WAVE = °/„d\n" , STANDING.WAVE) ; */ 
if ( ! STANDING.WAVE) 

{ getline(fp, inputline, 100); 

sscanf (inputline , "°/„lf", &tilt_angle) ; 
printf ("tilt_angle = %f \n" , tilt.angle) ; */ 

} 

else tilt_angle=0 . ; 

if ( ( ! STANDING_WAVE)&&(f abs (tilt.angle) >1 . e-10) ) 
nonzero_tilt=l ; 

else 

nonzero_tilt=0 ; 

getline(fp, inputline, 100); 

sscanf (inputline , "°/„d" , &num_detuning) ; 

for (i=0; i<num_detuning; ++i) 

{ getline(fp, inputline, 100); 

sscanf (inputline , "°/ lf", fedetuning [i] ) ; 

} 

getline(fp, inputline, 100); 

sscanf (inputline, "°/ lf", &INTERACTION_LENGTH) ; 

if (GAUSSIAN_X==0) INTERACTION.LENGTH =1.772; /* Sqrt(pi) 

getline(fp, inputline, 100); 

sscanf (inputline, "°/„d" , &CIRCULAR.BEAM) ; 

if (CIRCULAR.BEAM) 
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{ getline(fp, inputline, 100); 

sscanf (input line, "°/ lf", &D) ; 

} 

else 

{ getline(fp, inputline, 100); 

sscanf (inputline, '"/.If", &REC_HALF_HEIGHT) ; 

} 

getline(fp, inputline, 100); 

sscanf (inputline , "°/ lf", &ef f ective_V_to_apparent_V) ; 

if (GAUSSIAN_X==0) ef f ective_V_to_apparent_V = 1; 

/* For top hat, interaction region is equivalent to mode */ 

getline(fp, inputline, 100); 

sscanf (inputline , "°/ d" , &num_N_ef f ) ; 

for (i=0; i<num_N_eff; ++i) 

{ getline(fp, inputline, 100); 

sscanf (inputline, '"/.If", &N_ef f [i] ) ; 

} 

getline(fp, inputline, 100); 
sscanf (inputline , "°/ d" , &num_u_th) ; 

for (i=0; i<num_u_th; ++i) 

{ getline(fp, inputline, 100); 

sscanf (inputline , "%lf", &u_th[i]); 

} 

getline(fp, inputline, 100); 

sscanf (inputline, "°/ d M , &M0N0_VEL0CITY) ; 

if ( ! (M0N0_VEL0CITY) ) 

{ getline(fp, inputline, 100); 

sscanf (inputline, "°/„d", &USE_MAXWELL_B0LTZMANN) ; 
if ( ! (USE_MAXWELL_B0LTZMANN) ) 

{ getline(fp, inputline, 100); 

sscanf (inputline , "°/ lf", &_delta_v) ; 
getline(fp, inputline, 100); 
sscanf (inputline, "°/ lf", &_v0) ; 
} 

} 

getline(fp, inputline, 100); 

sscanf (inputline, "°/„d" , &GENERATE_PH0T0N_FILE) ; 
/* 

getline(fp, inputline, 100); 
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sscanf (input line , "°/ lf %lf", &a_coeff, &b_coeff); 
*/ 



getline(fp, inputline, 100); 

sscanf (input line, "°/„d" , &PUMPING.SIMULATED) ; 

if (PUMPING.SIMULATED) 

{ getline(fp, inputline, 100); 

sscanf (inputline, '"/.If °/„lf °/„lf", &PUMPING.LENGTH, 
&PUMP.WAIST, &PUMP.P0S); 
} 

getline(fp, inputline, 100); 
/* printf ("°/.s\n" , inputline) ; */ 
sscanf (inputline , "°/ s", summaryf ilename) ; 
if ( ! is alpha (summaryf ilename [0] ) ) 

summaryf ilename = SUMMARY_FILENAME_DEFAULT ; 
printf ("summaryf ilename = ,0 / s'\n", summaryf ilename) ; 

getline(fp, inputline, 100); 
printf (" : °/„s : \n" , inputline) ; 
if (isdigit (inputline [0] ) ) 
{ sscanf (inputline , "°/ d" , &n_seed) ; 

printf ("n_seed = °/ d\n" , n_seed) ; 

} 

getline(fp, inputline, 100); 
printf (" : %s : \n" , inputline) ; 
if (isdigit (inputline [0] ) ) 

{ 

sscanf (inputline, " / d" , &WAIT_FRACTION_TENTHS) ; 

printf ( "WAIT_FRACTION_TENTHS=°/ d\n" , WAIT_FRACTION_TENTHS) ; 

WAIT.FRACTION = (double) (( (double) WAIT_FRACTION_TENTHS) / 10.); 

/* WAIT.FRACTION = 5./10. ; */ 

printf ( "WAIT_FRACTION=°/„g\n" , WAIT.FRACTION) ; 

} 

getline(fp, inputline, 100); 
if (isdigit (inputline [0] ) ) 
{ sscanf (inputline, "°/„d" , &GAUSSIAN.X) ; 

printf ( "GAUSSIAN_X=°/.d\n" , GAUSSIAN.X) ; 

if (GAUSSIAN_X==0) INTERACTION.LENGTH =1.772; /* Sqrt(pi) */ 
if (GAUSSIAN_X==0) ef f ective_V_to_apparent_V = 1; 

> 
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/* 



getline(fp, inputline, 100); 
if (isdigit (inputline [0] ) ) 

sscanf (inputline, "°/„d" , &VARIABLE.FINESSE) ; 
printf ( "VARIABLE_FINESSE=°/ d\n" , VARIABLE.FINESSE) ; 



getline(fp, inputline, 100); 
if (isdigit (inputline [0] ) ) 

sscanf (inputline, "°/„d M , &ADJUST_N_EFF) ; 
printf ( "ADJUST_N_EFF=%d\n" , ADJUST_N_EFF) ; 



/* getline(fp, inputline, 100); */ 

/* printf ("This line should be blank: °/ s\n" , inputline); */ 
f close(fp) ; 

/* End of input from inputfile[K] */ 

/* compose output file names */ 

if ( (f p=f open(INDEX_FILE, "r"))==NULL) 
f ile_id=0; 
else 

{ fscanf(fp, "°/od", &f ile_id) ; 
f close(fp) ; 
} 

f ile_id++; 

f p=f open ( INDEX_FILE , "w"); 
fprintf(fp, "°/„d\n", f ile.id) ; 
f close (f p) ; 

sprintf (output.name, "°/„s.y„d", 0UTPUT.FILE, f ile.id) ; 



*/ 



sprintf (outputf ilename , "°/ s.out", inputf ile [K] ) ; 
printf ("outputf ilename = 'Zs'V, outputf ilename) ; 

sprintf (pnf ilename , "°/ s.pn", inputf ile [K] ) ; 
printf ("pnf ilename = '°/ s'\n", pnf ilename) ; 



/* 
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if (GENERATE_PHOTON_FILE) 
sprintf (photon_file_name, " /.s.°/„d", PHOTON_FILE, f ile.id) ; 

*/ 



if (GENERATE_PHOTON_FILE) 
sprintf (photonf ilename , "%s.pho", inputf ile [K] ) ; 

sprintf (statusf ilename , "°/„s.stat", inputf ile [K] ) ; 

time (&tp) ; 

sprintf (output_time , "°/„s", ctime(fetp)) ; 

f or (g_multiplier=0 . 1 ; g_multiplier <= 4.0; g_multiplier += 0.1) 

/* f or (g_multiplier=2 . 27; g_multiplier <= 2.3; g_multiplier += 0.1) */ 

{ printf ("g_multiplier = °/„g\n" , g_multiplier) ; 



/* calculate real world parameters */ 

zO=sqrt(rO*L/2. ) ; /* in cm */ 
fsr=3.e4/2./L; /* in MHz */ 

_w0=sqrt(lambda*10.*z0/PI) ; /* in micron */ 

f=3./4.* S qr(lambda/PI/_w0)*l .e-6; /* finesse */ 

_Gc_over_2pi=f sr/f inesse/1 . e6 ; /* in MHz */ 



_g_over_2pi=sqrt(2./PI*f sr*f *Ga_over_2pi/1000. ) ; /* in MHz */ 
_beam_radius=D/2./_wO; /* in _w0 */ 

if (! (STANDIMG.WAVE)) 
{ 

g_over_2pi /= 2; /* traveling wave correction - changed from */ 

} /* sqrt(2) to 2 by CFY */ 



_g_over_2pi *= g_multiplier ; 



/* parameters normalized to "2g" */ 
_Ga_over_2g=Ga_over_2pi/1000 . /2 . /_g_over_2pi ; 

/* find fastest rate and determine dtO */ 
X=(_g_over_2pi > _Ga_over_2g) ? _g_over_2pi : _Ga_over_2g; 
X=(X > _Gc_over_2pi) ? X : _Gc_over_2pi ; 
dt0=_g_over_2pi/X*dt0_init ; 

printf ("**** X=°/„g, _g_over_2pi=°/ g, dt0_init=°/„g\n" , X, 
_g_over_2pi, dt0_init) ; 



/* #ifndef LINUX */ /* amj */ 
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/* ef f ective_V_to_apparent_V= 1; */ 

/* amj change (double) volume_ratio() to 1 jul 14,98 */ 

/* CFY comments it out */ 
/* used in delta_t and <N> calculations */ 
/* #endif */ /* amj */ 

/* generate header information */ 
for (i=0; i<3; ++i) 
{ if (i==0) 

{ f p=f open(outputf ilename , "w") ; 
fprintf (fp, "# This file was generated by qt°/„d on °/„s\n" , 

VERS I ON , output _t ime ) ; 

} 

else 

if (i==l) 

{ fp=fopen(QT_LOG, "a"); 
fprintf (fp, 

"#Qt / d-°/ o d generated °/„s with the following parameters on %s\n#\n" , 
VERSION, REVISION, outputf ilename , output _t ime ) ; 

} 

else 

if (i==2) 

f p=stderr ; 

fprintf (f p , " ###########################################################\ n 11 ) ; 
fprintf (fp, 

"# Ga/2pi=%.llf kHz lambda=% . llf nm rO=°/„.llf cm\n" , 
Ga_over_2pi, lambda, rO) ; 

fprintf (fp, "# length=°/ . 31f cm f inesse=°/„ . llf million\n" , L, finesse); 
fprintf (fp, "# zO=%.31f cm wO=°/ .llf micron f =% . 21e\n" , 
zO, _w0, f); 

fprintf (fp, "# Gc/2pi=°/ .61f MHz g/2pi=°/„.31f MHz\n" , 
_Gc_over_2pi , _g_over_2pi) ; 

fprintf (fp, "# max_photon_size=°/ d max_atoms =0 / d\n" , 
_max_photon_size , _max_atoms) ; 

fprintf (fp, "# 2g*dt (initial) =%. 21e tot_atom_count=%d\n" , 
dtO, t ot _ at om_ count) ; 

fprintf (fp, "# eff ective_V_to_apparent_V=°/„.31e\n" , ef f ective_V_to_apparent_V) ; 

fprintf (fp, "# Gaussian along x & y\n"); 
if (STANDING. WAVE) 

fprintf (fp, "# stang-wave mode along z\n"); 

else 

fprintf (fp, "# traveling-wave mode along z with tilt=°/„.llf mrad 
(nonzero_tilt=°/ d) . \n" , tilt_angle, nonzero_tilt) ; 
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fprintf (fp, "# interaction length along x=%.31f\n", INTERACTION.LENGTH) ; 
if (CIRCULAR.BEAM) 

fprintf (fp, "# Circular atomic beam w/ diameter=°/„.31f \n" , D) ; 
else 

fprintf (fp, "# Rectangular atomic beam w/ beam half height=°/> . 31f \n" , 
REC_HALF_HEIGHT) ; 

if ( ! (MONO.VELOCITY) ) 
{ 

if (USE_MAXWELL_BOLTZMANN) 
fprintf (fp, 

"# Maxwell-Boltzmann velocity distr is used.\n"); 
else 

fprintf (fp, 

"# Lorentzian convoluted with Maxwell-Boltzman velocity distr is used 

with dv(FW)=°/ .31f and velocity peak v0=% . 31f \n" , _delta_v, _v0) ; 

> 

else fprintf (fp, "# Mono velocity is assumed. \n") ; 



#ifdef REUSE_C_ARRAY 

fprintf (fp, "# reuse C[] array=°/ d\n#\n" , USE_LAST_C_ARRAY) ; 
#endif 

if ( PUMP I NG_ S I MUL ATED ) 
{ 

/* _omega_Rabi=dimensionless , in 2g unit, pi pulse */ 

_omega_Rabi=PI/_t_int/PUMP_WAIST; /* CFY un-comments this line */ 
if ( (INTERACTI0N_LENGTH+PUMPING_LENGTH)/2 . > PUMP.POS) 
PUMPING_EFFECT_0N=1 ; 
else PUMPING_EFFECT_ON=0 ; 

PUMPING_EFFECT_ON=0 ; /* ******** CFY for test ****** */ 
fprintf (fp, "# Pumping process is simulated: \n") ; 

fprintf (fp, "# pumping length=°/ . 31f pump waist=°/„.31f pump pos=% . 31f \n" , 
PUMPING.LENGTH, PUMP.WAIST, PUMP.POS); 

> 

if ( INCLUDE_ATOM_DECAY) 

fprintf (fp, "# Atom spontaneous emission included. \n") ; 
if ( INCLUDE_D_STATE_DECAY) 
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fprintf(fp, "# D-state branching included. \n") ; 

#ifdef USE_RK4_ALG0RITHM 

fprintf(fp, "# RK4 algorithm used.\n"); 

#endif 

if (i<2) fclose(fp); 
} 



if (GENERATE_PHOTON_FILE) 
{ 

fp=fopen(photonfilename, "w"); 

fprintf(fp, "# This file was generated by qt°/„d on °/„s\n" , 
VERSION, output _t ime ) ; 
f close (f p) ; 
} 

/* Check array size */ 
max_k=two_to (_max_atoms) ; 
C_size=(long) (_max_photon_size) *max_k; 

total_array_size=C_size*sizeof (double) *8 /* display purpose */ 
+(long) ( (_max_atoms+l)*A_INFO_ITEMS)*sizeof (double) 
+(long) (_max_atoms+l) *max_k*sizeof (double) ; 
if (total_array_size > MAX.MEMORY) 

{ fprintf (stderr , "# Too large array size (%ld kB) . 
Reduce it below '/.Id kB.\n\ total_array_size/1024L, MAX_MEM0RY/1024L) ; 
exit(l) ; 
} 

/* reserve memory space for _C[] */ 

_C=(double *)malloc (C_size*sizeof (double) ) ; /* real part only */ 

_D=(double *)malloc(C_size*sizeof (double) ) ; /* imaginary part only */ 

_C_l=(double *)malloc(C_size*sizeof (double) ) ; 

_D_l=(double *)malloc(C_size*sizeof (double) ) ; 

_Pn_dist= (double *)malloc (_max_photon_size*sizeof (double) ) ; 

_Pn_dist_mean= (double *)malloc (_max_photon_size*sizeof (double) ) ; 

_Pn_dist_mean_beginning= (double *)malloc (_max_photon_size*sizeof (double) ) ; 

_detuning_sum= (double *)malloc (max_k*sizeof (double) ) ; 

_bits= (double **)malloc(max_k*sizeof (double) ) ; 

for (i=0; i<max_k; ++i) 

_bits [i] = (double *)malloc (_max_atoms*sizeof (double) ) ; 

#ifdef USE_RK4_ALG0RITHM 

_C_2=(double *)malloc(C_size*sizeof (double) ) ; 
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_C_3=(double *)malloc(C_size*sizeof (double) ) ; 
_D_2=(double *)malloc(C_size*sizeof (double) ) ; 
_D_3=(double *)malloc(C_size*sizeof (double)) ; 
#endif 

for (i=0; i<2; ++i) 

{ if (i==0) fp=f open(outputf ilename, "a"); 
else fp=stderr; 

fprintf(fp, "# total array size=%ld kBytes\n" , 



/* initialize rand() */ 

#ifdef RANDOM_INIT 
srand( (int)time (&tp) ) ; 

#else 
srand(l) ; 

#endif 

#ifdef REUSE_C_ARRAY 
if (USE_LAST_C_ARRAY) read_status_f ile () 
#endif 



for (m=0; m<num_N_ef f ; ++m) { N_eff [m] loop *********/ 

for (i=0; i<3; ++i) 

{ if (i==0) fp=fopen(QT_LOG, "a"); 
else if (i==l) f p=f open(outputf ilename , "a"); 
else if (i==2) fp=stderr; 
f printf (f p, 



total_array_size/1024L) ; 



if (i==0) fclose(fp); 
} 



"# 

f printf (f p, 



\n"); 



# %s\n\ inputfile[K]) ; 



fprintf(fp, "# <N>=°/ .31f \n", N_ef f [m] ) ; 
if (i<2) fclose(fp); 



} 



for (mm=0; mm<num_u_th ; ++mm) 
{ _u_th=u_th [mm] ; 

if (USE_MAXWELL_BOLTZMANN) 
_vO=l; 



****** u_th [mm] loop *********/ 



/* amj june 15, 98 */ 
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if (MONO.VELOCITY) 
_vO=l; 



/* amj jul 14, 98 */ 



if (GAUSSIAN.X) 
t_int=sqrt(PI)*_wO/(_vO*_u_th) ; 

/* Based on Rabi oscillation, not on mode volume (smaller by sqrt(2)). 
/* In microsec */ 

else t_int= sqrt(PI)*_wO/(_vO*_u_th) ; 
/* amj change 2 . *_wO/(_vO*_u_th) ; to sqrt(PI)*_wO/(_vO*_u_th) ; jul 14,98 
/* amj multiply _u_th by _v0 June 15, 98 _v0 is peak of velocity */ 
/* distribution in units of _u_th */ 

/* parameters normalized to "2g" */ 
_t_int=t_int*2 . *_g_over_2pi*2 . *PI ; 

if (VARIABLE.FINESSE) 

{ 

_Gc_over_2pi=(_u_th/u_th[0] ) * f sr/f inesse/1 . e6; /* in MHz */ 
printf ("**** _Gc_over_2pi=y g\n" ,_Gc_over_2pi) ; 

} 

/* Doppler shift due to tilt for v=u_th in 2g unit, including zero tilt */ 

mean_shift=(_v0*_u_th)/lambda*sin(tilt_angle/1000.)*l000. ; 

/* in MHz */ 

_mean_shif t_over_2g=mean_shif t/ (2 . *_g_over_2pi) ; 
/* amj change _u_th to _vO*_u_th June 15, 98 */ 

/* _v0 is peak of velocity distribution */ 

for (i=0; i<3; ++i) 

{ if (i==0) fp=fopen(QT_LOG, "a"); 
else if (i==l) f p=f open(outputf ilename , "a"); 
else if (i==2) fp=stderr; 
f printf (f p, 

m # \ n „) 

fprintf(fp, "# _u_th=°/ . 31f t_int=%.31f microsec\n" , 

_u_th, t_int) ; 

if (nonzero_tilt) 

fprintf(fp, "# mean transverse Doppler shif t=°/ . 31f \n" , 
mean_shift) ; 

if (i<2) fclose(fp); 

} 
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/* delta_t is the mean time interval between atoms in microsec */ 
/* based on weighted volume integral */ 



delta_t=t_int/N_ef f [m] ; 
/* amj change */ 

/* INTERACTION_LENGTH*eff ective_V_to_apparent_V*_wO/(_vO*_u_th)/N_ef f [m] ; 

/* to t_int/N_ef f [m] ; jul 14,98 */ 

/* amj change _u_th to _vO*_u_th June 15, 98 */ 

/* _v0 is peak of velocity distribution in units of _v_th*/ 

/* normalized to 2g */ 
_Delta_t=delta_t*2 . *PI*_g_over_2pi*2 . ; 

for (i=0; i<2; ++i) 

{ if (i==0) f p=f open(outputf ilename , "a"); 
else fp=stderr; 
f printf (f p, 

"# 2g*t_int=°/ .21e 2g*Delta_t=°/ . 21e\n" , _t_int, _Delta_t) ; 
if (i==0) fclose(fp); 



for (mmm=0; mmm<num_detuning ; ++mmm) 
{ /** detuning [mmm] loop **/ 

/* readjust dt considering detuning effect */ 

if (nonzero_tilt) 

{ /* detuning [] measured in _mean_shif t_over_2g unit */ 
_CA_detuning=detuning [mmm] *_mean_shif t_over_2g ; 
dt=dtO/sqrt (sqr (_CA_detuning-_mean_shif t_over_2g)+ 1.) ; 
printf ("dt=%g\n\dt) ; 
} 

else 

{ /* detuning [] measured in MHz */ 
_CA_detuning=detuning [mmm] /(2*_g_over_2pi) ; 
dt=dtO/sqrt (sqr (_CA_detuning) +1 . ) ; 



for (i=0; i<3; ++i) 

{ if (i==0) fp=fopen(QT_LOG, "a"); 
else if (i==l) fp=fopen(outputf ilename, "a"); 
else if (i==2) fp=stderr; 
f printf (fp, 

M # \ n „) . 

if (nonzero_tilt) 
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fprintf (fp, 

"# detuning=%g (v/lambda*theta) or °/ g MHz dt (initial) =°/„g\n\n" , 
detuning [mmm] , detuning [mmm] *mean_shift , dt) ; 
else 

fprintf (fp, "# detuning=°/ g MHz dt (initial) =°/ g\n\n" , 
detuning [mmm] , dt) ; 
fprintf (fp, 

"# t n_aver n2_aver Mandel_Q P_N_mode atom_count\n") ; 

/* amj change dn_aver to Mandel_Q and add n2_aver */ 

if (i<2) fclose(fp); 

} 

/* start benchmark clock */ 
Tl=clock() ; 

/* initialize */ 

for (i=0; i<=_max_atoms ; ++i) P_atom[i]=0; 

for (i=0; i<_max_photon_size ; ++i) _Pn_dist_mean[i]=0. ; 

for (i=0; i<_max_photon_size ; ++i) _Pn_dist_mean_beginning[i]=0. ; 

if ( ! (USE_LAST_C_ARRAY) ) { 

/* Initialize _C[] arrays to zero */ 

for (ii=0; ii<C_size; ++ii) {_C[ii]=0.; _D[ii]=0.;} 

/* initially vacuum + no atoms */ 
/* _C[0]=1. ; */ 

_C[n_seed] =1.; /* Seed with n_seed photons! */ 

/* initialize _C[] using one-atom micromaser theory 

C_f rom_micromaser_theory (N_ef f [m] , a_coeff, b_coeff); */ 

_M_atom=0; 

t=0.; /* simulation starts from t=0 by default. If last C[] used, 

t is the last time in the previous simulation. */ 

} 

t=0.; /* initial time of new simulation */ 
Ti=t; 

next_entry_time=t ; 

/* time averaging begins when t>Ti2 */ 

Ti2=Ti+_Delta_t* (double) (tot_atom_count) *WAIT_PORTION_BEFORE_AVERAGING ; 
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Ti2new = Ti+_Delta_t* (double) (tot_atom_count) * 0.1; 

/* initialize count variables */ 
atom_count=0; /* number of atoms injected */ 

atom_count_since_Ti2=0; /* number of atoms injected since t=Ti2 */ 
photon_count=0; /* number of photons emitted */ 
too_many_atom_count=0 ; /* number of atoms fail to be injected */ 
num_excited_injected_atoms=0 . ; /* sum of excited state probabilities 

of injected atoms */ 
num_excited_exiting=0 ; /* number of exiting atoms in excited state */ 

if ( INCLUDE_ATOM_DECAY) spont_count=0 ; /* number of spont emission decays */ 
if (INCLUDE_D_STATE_DECAY) d_state_decay_count=0 ; 

/* initialize averaged variables */ 

n_mean=0.; /* amj This is supposed to initialized */ 

n2_mean=0 . ; /* amj add */ 

Mandel_Q=0.; /* amj add */ 

n_aver=0.; /* time averaged mean photon number */ 
n2_aver=0.; /* amj change dn_ave to n2_ave */ 

/* time averaged variance of photon number */ 

n_sum=0 . ; 

n2_sum=0.; /* amj change dn_sum to n2_sum */ 

emission_prob=0 . ; /* averaged ratio number of emitted photons 
to total number of injeced atoms */ 

vacuum_time=0 . ; 
num_dt_steps=0 ; 

/* simulation starts at t=Ti, averaging begins at t=Ti2 */ 
do { 

if (t>=next_entry_time) { 

excited_state_prob=inject_new_atom(t) ; /* amj */ 
if (excited_state_prob<0 . ) 

{ too_many_atom_count++; 

printf("** count = °/„d out of °/„d (°/ .2g / °/ )\n" , too_many_atom_count , atom_ 

} 

/* inject_new_atom() returns negative value if not successful to inject an atom. 

Otherwise, returns excited state probability */ 
else if (excited_state_prob>99) atom_decay_to_d_count++; /*amj if decay to d state 
else { /* successful injection */ 

atom_count++; /* number of successfully injected atoms */ 
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/* fprintf (stderr, "."); amj */ 
if (t>=Ti2) { 

atom_count_since_Ti2++ ; 

num_excited_inj ected_atoms+=excited_state_prob ; 

} 

/* user can stop execution */ 
#ifdef CLICK. INTERRUPT /* Macintosh specific */ 
if (key_pressed() ) { 

fprintf (stderr , "Stop now? (stop=l, continue=0)\n") ; 
C scanf ("°/ d" , festop.now) ; 

if (stop_now) break; 

} 

#endif 
} 

_max_k=two_to (_N_atom) ; /* update _max_k */ 
next_entry_time = t + get_arrival_time() ; 

if ( DEBUG_LEVEL> 1 ) 

printf ("N=°/„d created at t=°/„lf and will exit at t=°/,lf . 
Next new arrival at t=°/„lf \n\n" , 

_N_atom, t, _A_inf o [_N_atom] [2] , next_entry_time) ; 
} /* end of adding a new atom */ 

if (_N_atom>0) { 

/* coherent evolution of wavefunction */ 
evolve_wave_f unction(dt , t) ; 
/* only when there is at least one atom */ 
t+=dt ; 

/* take care of atomic spontaneous emission */ 

if (INCLUDE.ATOMJDECAY) { 

if (t<Ti2) how_many_atoms_decay (dt) ; 

else spont_count+=how_many_atoms_decay(dt) ; 

} 

/* check for D- state decay */ 

if (INCLUDE_D_STATE_DECAY) { 

if (t<Ti2) how_many_d_state_decay (dt) ; 

else d_state_decay_count+=how_many_d_state_decay (dt) ; 

/* once decayed to D state the 

atom is not interacting with field */ 

> 

/* check for exiting atom */ 
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for (i=l; i<=_N_atom;) { 

if (_A_inf o [i] [2] < t) { /* remove ith atom */ 
if ( DEBUG_LEVEL> 1 ) 

printf ("Removing °/ d-th atom at t=%lf out of °/ d\n" , i, t, _N_atom) ; 

if (t<Ti2) remove_atom(i) ; 

else num_excited_exiting+=remove_atom(i) ; 

/* remove_atom() returns 1 if 

the atom leaves the cavity in 
the excited state */ 

_max_k=two_to(_N_atom) ; /* update _max_k */ 
> 

else i++; 
> 

} /* if _N_atom>0 */ 

else t+=dt; /* no evolution if no atom there */ 
num_dt_steps++ ; 

/* check for cavity decay */ 

if ((_N_atom==0)&&(n_mean<VERY_SMALL)) { /* complete vacuum */ 
big_dt=next_entry_time-t ; 
vacuum_time+=big_dt ; 
P_atom[0]+=(long) (big_dt/ (float) dt) ; 

/*amj put float to avoid getting smaller values. CFY: ??? */ 
_Pn_dist_mean [0] +=big_dt/dt ; 
big_dt=big_dt/ (float) dt; 

t+=(long)big_dt*dt ; /* amj */ /* nothing happens until 
next atom comes in */ 
num_dt_steps+=(long)big_dt ; /* amj */ 
} ' 

else { /* not vacuum, equal stepsize dt */ 

if (cavity_decay_occurs (dt , &n_mean, &n2_mean)) 
{/* amj add &n2_mean */ 

if (t>=Ti2) photon_count++; 
if (GENERATE_PH0T0N_FILE && 
(atom.count > tot_atom_count * WAIT.FRACTION )) { 
f p=f open(photonf ilename , "a"); 
fprintf(fp, "°/ d °/„lf \n" , photon_count , t) ; 
f close (f p) ; 

} 

} 

/* if (atom_count == (floor) 
((f loat)tot_atom_count * WAIT.FRACTION) && t>=next_entry_time) */ 
n_sum+=n_mean ; 
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/* printf ("n_mean = °/„g\n" , n_mean) ; 
printf ("n_sum = °/ g\n" , n_sum) ; */ 

n2_sum+=n2_mean; /* amj using average over "10" atoms */ 

/* P_atom is updated from t=Ti, not Ti2 */ 
P_atom[_N_atom] ++; /* meaningful for equal size steps */ 

} 



if (t>=Ti2) { 

/* cumulative averag */ 

/* n_aver=n_sum/(t-Ti2) ; amj */ 

/*n2_aver=n2_sum/(t-Ti2) ;*/ /* amj change dn_aver to n2_aver 
and d2_sum to n2_sum */ 
emission_prob= (double) photon_count/ (double) atom_count_since_Ti2; 
/* Mandel_Q=n2_aver/n_aver-n_aver-l . ; amj */ /* amj */ 

} 

/* if ((atom_count°/ 100 == && (atom_count<1000) ) ) 

printf ("atom count=°/ d\n" , atom_count) ; */ 

if ( (atom_count == 100) && (t>=next_entry_time) ) 
{ T2=clock(); 

time.elapsed = (T2-T1)/CL0CKS_PER_SEC; 

printf ("100 atoms in °/„g sec. Estimated time for °/„d atoms: 
°/ g sec (%.2g hours)\n", time.elapsed, t ot_ at om_ count , 
time_elapsed* (double) t ot _ at om_ count/ 100 , 
time.elapsed* (double)tot_atom_count/100/3600) ; 
} 

if (atom.count == (floor) ( (float) tot_atom_count * WAIT.FRACTION) 
&& t>=next_entry_time) 

{ printf ("Skipping the first °/„d atoms for P_n distribution and Q.\n", 
(int) (floor) ( (f loat)tot_atom_count * WAIT_FRACTI0N) ) ; 
for (i=0; i<_max_photon_size ; ++i) 

{ _Pn_dist_mean_beginning [i] = _Pn_dist_mean [i] ; 
/* printf ("%d °/„g\n", i, 

_Pn_dist_mean_beginning[i] /( (t-Ti)/dt) ) ; */ 
/* printf ("'/.g °/og\n", t,Ti); */ 
> 

} 

if (ADJUST_N_EFF && (atom_count == (floor) ( (float) tot_atom_ count 
* WAIT_FRACTI0N/2.) && t>=next_entry_time) ) 

{ printf ("Skipping the first %d atoms for calibrating N_atom.\n", 
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(int) (floor) ( (float) tot_ atom. count * WAIT_FRACTI0N/2 . ) ) ; 

/* Calculate intracavity atom number distribution */ 
P_sum=0.; /* sum of P_atom[] */ 
P_N_mean=0 . ; 

for (i=0; i<=_max_atoms ; ++i) { 
P_sum+= (double) P_atom[i] ; 
/* P_atom[] updated during above simulation */ 
> 

for (i=0; i<=_max_atoms ; ++i) 

P_N_mean+= (double) i*(double)P_atom[i] ; 
P_N_mean/=P_sum ; 

printf ("P_N_mode=°/,g\n" ,P_N_mean * ef f ective_V_to_apparent_V) ; 



if (ADJUST_N_EFF) 

{ new_N_eff = N_ef f [m] / 

(P_N_mean * ef f ective_V_to_apparent_V / N_eff[m]); 
delta_t=t_int/new_N_ef f ; 
_Delta_t=delta_t*2 . *PI*_g_over_2pi*2 . ; 
printf ("new_N_eff=°/og, delta_t=°/„g\n" , new_N_eff, delta.t) ; 

} 

if (AD JUST.N.EFF) 
{ 

_Gc_over_2pi=(_u_th/u_th[0] ) * (f sr/f inesse/1 . e6) * 

(P_N_mean * ef f ective_V_to_apparent_V / N_eff[m]); 
printf ("Gc multiplied by °/„g\n" , 
P_N_mean * ef f ective_V_to_apparent_V / N_ef f [m] ) ; 

printf ("**** _Gc_over_2pi=%g\n" ,_Gc_over_2pi) ; 

} 



#ifdef SHOW.PROGRESS 

if ( (DEBUG_LEVEL>0) I I ( ( ! (atom_count°/„1000) )&& 
(t>=next_entry_time) ) ) 

{ 

num_dt_aver_steps= (double) (num_dt_steps-num_dt_aver_steps_last) ; 
/* amj */ 

/* amj mun_dt_aver_steps # of dt over which n_aver is found */ 
n_aver=n_sum/num_dt_aver_steps ; /* amj */ 

n2_aver=n2_sum/num_dt_aver_steps ; /* amj */ 
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Mandel_Q=n2_aver/n_aver-n_aver-l . ; /* amj */ 

/* Calculate intracavity atom number distribution */ 
P_sum=0.; /* sum of P_atom[] */ 
P_N_mean=0 . ; 

for (i=0; i<=_max_atoms ; ++i) 
P_sum+= (double) P_atom [i] ; 

for (i=0; i<=_max_atoms ; ++i) 

P_N_mean+= (double) i* (double) P_atom [i] /P_sum ; 

for (i=0; i<2; ++i) 
{ 

if (i==0) f p=f open(outputf ilename , "a"); 
else fp=stderr; 

fprintfCfp, " °/„.61e °/..61e / .61e °/„+.61e °/„g °/„14d\n" , 
t, n_aver, n2_aver, Mandel_Q, 

ef f ective_V_to_apparent_V*P_N_mean, atom_count) ; 
/* output */ 
if (i==0) f close (fp); 

} 

/* printf ("num_dt_aver_steps=°/ g\n" , num_dt_aver_steps) ; 

printf ( "n_sum=°/ g\n" , n_sum) ; 

printf ("n_aver=°/ g\n" , n_aver) ; */ 

num_dt_aver_steps_last= (double) num_dt_steps ; 
/* amj seting for next time */ 

n_sum=0; /* amj */ 

n2_sum=0; /* amj */ 

} 

#endif 

/*if (! (num_dt_steps°/ (long) (0.05/dt))) 

printf ('7.1e °/„le %le\n N , t, sqr(_C[l] )+sqr(_D[l] ) , sqr (_C[3] )+sqr(_D [3] ) ) ; */ 

/* amj test starts jul 27, 98 */ 
/* 

for (i=0; i<2; ++i) 
{ 

if (i==0) 

fp=f open (outputf ilename , "a"); 
else fp=stderr; 
fprintfCfp, "'/.If ", t) ; 
for (ii=0; ii<=5; ++ii) 
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fprintf(fp, "'/.If 7.1f ", _C[ii], _D[ii]); 
printf ("\n") ; 
if (i==0) f close (fp); 

} 

*/ 

/* amj test endss jul 27, 98 */ 



} while (atom_count<=tot_atom_count) ; 

/* ******* END OF SIMULATION RUN ********** */ 
f printf (stderr, "\n"); 

atom_count — ; /* simulation stops at last atom entry time */ 
Ti=0. ; 

if (DEBUG_LEVEL>1) { 

printf ("*** Contents of _A_inf o [] . . . \n") ; 
for (ii=l; ii<=_N_atom; ++ii) { 

for (i=0; i<6; ++i) printf ("'/.If " , _A_inf o [ii] [i] ) ; 

printf ("\n") ; 

} 

printf ("*** Contents of _C[]...\n"); 
for (i=0; i<=l; ++i) { 
for (ii=0; ii<=4; ++ii) 

printf ("%lf ", _C [ii*_max_photon_size+i] ) ; 

printf ("\n") ; 

} 

} 

/* Calculate intracavity atom number distribution */ 
P_sum=0.; /* sum of P_atom[] */ 
P_N_mean=0 . ; 

for (i=0; i<=_max_atoms ; ++i) { 
P_sum+= (double) P_atom [i] ; 

/* P_atom[] updated during above simulation */ 
} 

for (i=0; i<=_max_atoms ; ++i) P_N_mean+= (double) i* (double) P_atom [i] ; 
P_N_mean/=P_sum ; 

/* Subtract _Pn_dist_mean_beginning here */ 
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for (i=0; i<_max_photon_size; ++i) _Pn_dist_mean[i] -= 
_Pn_dist_mean_beginning [i] ; 

/* global variable _Pn_dist_mean [] was updated since t=Ti */ 
for (i=0; i<_max_photon_size ; ++i) 
_Pn_dist_mean [i] /= (t-Ti) * (1-WAIT_FRACTI0N) /dt ; 

/* stop benchmark clock */ 
T2=clock() ; 



/* averaging stuff taken from "SHOW.PROGRESS" routine */ 
num_dt_aver_steps= (double) num_dt_steps-num_dt_aver_steps ; 
/* amj num_dt_aver_steps # of dt over which n_aver is found */ 
n_aver=n_sum/num_dt_aver_steps ; 
n2_aver=n2_sum/num_dt_aver_steps ; 
Mandel_Q=n2_aver/n_aver-n_aver-l . ; 

/* Notifying key informations on data */ 

for (i=0; i<3; ++i) { 

if (i==0) fp=fopen(QT_LOG, "a"); 

else if (i==l) fp=f open(outputf ilename, "a"); 

else fp=stderr; 

fprintf(fp, "# °/„lf sec taken... \n\ (T2-T1)/CL0CKS_PER_SEC) ; 
fprintf(fp, "# °/„d atoms successfully injected and %d atoms 
used for time averaging. \n" , 
atom_count, atom_count_since_Ti2) ; 

fprintf(fp, "# °/„d atoms ( / g /o°/o) were neglected due to overflow 
of atoms during entire simulation\n" , 
too_many_atom_count , 

100* (float) too_many_atom_count/ (float) atom_count) ; 
fprintf(fp, "# °/„d photons decayed during 2gT=°/„g (averaging time)\n", 
photon_count , t-Ti2) ; 

fprintf(fp, "# %d atoms exit the cavity in excited state during 
entire simulation. \n" , 
num_excited_exiting) ; 

fprintf(fp, "# Effective number of excited atoms was °/ g during 
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entire simulation. \n" , 

num_excited_injected_atoms) ; 

fprintf(fp, "# Pumping efficiency was °/ g\n" , 

num_excited_inj ected_atoms/ 

(double) atom_count_since_Ti2) ; 

if ( INCLUDE_ATOM_DECAY) 

fprintf(fp, "# %d atoms undergone spontaneous emission during averaging\n" , 
spont_count) ; 

if ( INCLUDE_D_STATE_DECAY) 

fprintf (fp, "# °/„d atoms decayed to 1D2 state during averagin\n" , 
d_state_decay_count) ; 

fprintf (fp, "# <n>=°/„g Mandel Q=%g emission_prob=°/„g\n" , 
n_aver, Mandel_Q, emission_prob) ; 

Gc_f inal [mmm] = (double) photon_count/n_aver/ (t-Ti2)*2. *_g_over_2pi ; 
fprintf (fp, 

"# From simulation, 2g*Delta_t (sim)=%lf and Gc/2p(sim)=%lf MHz\n" , 
(t-Ti) / (double) atom_count , Gc_f inal [mmm] ) ; 
fprintf (fp, "# During simulation, the cavity was completely 
empty for °/„.31f percent of time.\n", vacuum_time/(t-Ti)*100. ) ; 

fprintf (fp, "# P(N) distribution averaged for entire simulation\n") ; 
fprintf (fp, "# Compared with Poisson distribution\n") ; 
p_poisson = exp(-l*N_ef f [m] /ef f ective_V_to_apparent_V ); 
for (j=0; j<=_max_atoms ; ++j) 
{ if(j!=0) 

p_poisson *= (N_ef f [m] /ef f ective_V_to_apparent_V) / (float) j; 
fprintf (fp, "°/ d °/„lf °/ lf\n", j, (double) P_atom [j] /P_sum, p_poisson) ; 

} 



/* fprintf (fp, "# P_n distribution averaged for simulation\n") ; */ 
final_n_avg = 0; 
final_n2_avg = 0; 
f inal_normalize = 0; 

for (j=0; j<_max_photon_size; ++j) 

{ /* fprintf (fp, "°/„d °/„lf \n" , j, (double) _Pn_dist_mean [j] ) ; */ 
f inal_n_avg += (double) j * (double) _Pn_dist_mean[j] ; 
final_n2_avg += (double) j * (double) j * (double) _Pn_dist_mean [j] ; 
f inal_normalize += (double) _Pn_dist_mean[j] ; 
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} 



final_n_avg /= f inal_normalize; 
final_n2_avg /= f inal_normalize ; 

fprintf (fp, "\n# <n> = °/„g, <n~2>=°/ g, Q=°/ g\n\n" , 

final_n_avg, f inal_n2_avg, (f inal_n2_avg - 

f inal_n_avg*f inal_n_avg - f inal_n_avg)/f inal_n_avg) ; 

fprintf (fp, "# normalization was Sum(P_n) = °/ g\n" , f inal_normalize) ; 

fprintf (fp, "# Averaged apparent number of atoms in 
BEAM_AREA* INTERACT! ON_LENGTH=°/ g\n " , P_N_mean) ; 

P_N_mode=ef f ective_V_to_apparent_V*P_N_mean; 

fprintf (fp, "# Averaged effective number of atoms in mode 
volume=°/ g\n\n" , P_N_mode) ; 

if (i<2) fclose(fp); 

n_f inal [mmm] =n_aver ; 
Q_f inal [mmm] =Mandel_Q ; 
N_f inal [mmm] =P_N_mode ; 
} 

/* SAVE SUMMARY INFORMATION */ 

fp = f open(summaryf ilename , "a"); 
time (&tp) ; 

sprintf (output_time , "°/„s", ctime(&tp)); 

fprintf (fp, "# °/„s °/ s", inputf ile [K] , output_time) ; 
fprintf (fp, 

"# N_ef f , Finesse, N_eff * Finesse, <n>, Q, det, u_th, N_ex, 
theta_pump, g_mult\n"); 

printf ("_Gc_over_2pi=°/ g, t_int=°/ g \n", _Gc_over_2pi , t_int) ; 
N_ex = P_N_mode / ( _Gc_over_2pi * 2*pi * t_int) ; 
theta_pump = sqrt(N_ex) * _t_int / 2.; 

fprintf (fp, "°/.g %g %g %g %g °/.g %g %g %g °/.g\n", 
P_N_mode, finesse, P_N_mode * finesse, 
final_n_avg, (f inal_n2_avg - f inal_n_avg*f inal_n_avg - 
f inal_n_avg)/f inal_n_avg, detuning [mmm] , 
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_u_th, N_ex, theta_pump, g_multiplier) ; 
f close (f p) ; 

/* WRITE PN DISTRIBUTION */ 

fp = f open(pnf ilename , "w"); 
for (j=0; j<_max_photon_size; ++j) 

{ fprintfCfp, "°/„d °/„lf \n" , j, (double) _Pn_dist_mean [j] ) ; 
} 

f close (f p) ; 

printf("P_n distribution written to ' / s'\n M , pnf ilename) ; 
} 

detuning [mmm] loop ******************/ 
/* generate summary sheet */ 
if (num_detuning>l) { 
for (i=0; i<3; ++i) { 
if (i==0) fp=fopen(QT_LOG, "a"); 
else if (i==l) fp=fopen(outputf ilename , "a"); 
else fp=stderr; 

fprintfCfp, "detuning <n> Q G_c <N>\n"); 

for (mmm=0; mmm<num_detuning ; ++mmm) { 

fprintfCfp, '7„.31e °/ .31e %.31e %.31e %.31e\n", 

detuning [mmm] , n_f inal [mmm] , Q_f inal [mmm] , Gc_f inal [mmm] , 

N_f inal [mmm] ) ; 

} 

if (i<2) fclose(fp); 

} 

} 

y z***************** N_u_th[mm] loop ******************/ 
y N_ef f [m] loop ******************/ 

/* write_status_f ile() ; */ 

/* 

free(_C) ; 

free(_D) ; 

f ree(_Pn_dist) ; 

f ree (_Pn_dist_mean) ; 

free (_bits) ; 
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f ree (_detuning_sum) ; 

CFY: free() is giving Seg Fault for higher N 

*/ 



/* 

#ifdef USE_RK4_ALG0RITHM 

free(_C_l) ; 

free(_C_2) ; 

free(_C_3) ; 

free(_D_l) ; 

free(_D_2) ; 

free(_D_3) ; 

#endif 

*/ 

} 

} /* **************** GRAND LOOP ****************** */ 



fprintf (stderr, "*** End of Run ***\n"); 

return 1 ; 

> 

void testroutine(void) 
{ 

int hist [100] ; 

int i, j, bin, npoints = 10000; 
double arrt ; 

_Delta_t =1.0; 

for(i=0;i<100;i++) 
hist[i]=0; 

f or ( j =0 ; j <npoints ; j ++) 

{ arrt = get_arrival_time() ; 
bin = f loor((arrt*20)) ; 
hist [bin] ++; 

/* printf ("hist[°/„d] =°/.d\n" , bin, hist [bin] ) ; */ 
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} 



for(bin=0;bin<100;bin++) 

{ arrt = (double) bin /20.; 
printf("°/.g °/„g °/„g\n", arrt, ( (double) (hist [bin] )) / 
(npoints*_Delta_t/20.) , exp(0 - arrt/_Delta_t)/_Delta_t ); 
} 

exit(O) ; 

} 

int getline(FILE *ff, char *s, int lim) 
{ int c , i ; 

for(i=0; i<lim-l && (c=getc (f f ) ) ! =E0F && c!='\n'; ++i) 

s[i] = c; 
if (c== '\n') 
{ s[i] = c; 
++i; 

} 

s[i] = '\0'; 
return i ; 

} 

double get_arrival_time() 
{ double gresult; 

gresult= _Delta_t * log(l. 0/(1.0 - drand())); 
return gresult; 

} 



#ifdef REUSE_C_ARRAY 
#include " status. c" 
#endif 

#include "newatom2-d. c" 
#include "exitatom2-d. c" 
#include "decay2-d.c" 
#include "evolve2-2-d. c" 
#include "spont2-d.c" 
#include "pumping-d. c" 
#include "d_state2-d. c" 
#ifdef CLICK. INTERRUPT 
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#include "key_pressed. c" 
#endif 

/* #ifndef LINUX */ /* amj */ 

#include "volume . c " 

/* #endif */ /* amj */ 



B.2 newatoms2-d.c 

#define DELTA_T_MAX 10. /* in _Delta_t scale */ 
#define V_MIN 0.0 
#define V_MAX 4. 
#define A 2.43952 

double get_velocity() ; 

float get_height() 
{ 

double a, b, bb, y; 

if ( ! (CIRCULAR_BEAM) ) y=REC_HALF_HEIGHT*drand() ; 
if (CIRCULAR.BEAM) { 
do { 

a=drand() ; 

b=sqrt (1 . -sqr (a) ) ; 

bb=drand() ; 

} while (bb>b) ; 

y=_beam_radius*a; 

> 

return y; 
> 

/* inject_new_atom() returns negative value if not successful to inject an atom 
Otherwise, returns excited state probability */ 

double inject_new_atom (double t) 
{ 

double v, y, excited_prob; 
long m, max_m; 

void construct_C() ; 
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long two_to(int); 

/* amj #ifdef PUMPING. SIMULATED */ 
void get_state_amplitude (double) ; 
/* amj #endif */ 

_N_atom++; 

/* if(_N_atom == 1) printf ("First atom enters at time °/„g\n", t) ; */ 

/* printf ("'/.d ", _N_atom) ; */ 

if (_N_atom > _max_ atoms) { 

f printf (stderr, "Too many atoms! "); 

fprintf (stderr, "t=°/„lf ", t) ; 

_N_atom — ; /* throw away this new atom */ 

excited_prob = -1.; /* negative means failure */ 

goto LAST.LINE; 

} 

/* get new atomic parameters */ 
if (MONO.VELOCITY) v=l.; 
else v=get_velocity () ; 

/* amj #ifdef PUMPING.SIMULATED *//* amj move this block to this position */ 

if (PUMPING.SIMULATED) /* amj */ 
{ 

/* check overlap, include spont . emiss. 

PUMP.POS, INTERACTION_LENGTH , PUMPING.LENGTH to be used */ 

/* _omega_Rabi=dimensionless , in 2g unit, pi pulse */ 
_omega_Rabi=PI*v/_t_int/PUMP_WAIST; /* amj */ 

if ( (PUMP_P0S+PUMPING_LENGTH/2 . ) <= INTERACTI0N_LENGTH/2 . ) 
{ 

_A_inf o [_N_atom] [4] =1 . ; /* ground state */ 
_A_inf o [_N_atom] [5] =0 . ; /* excited state */ 
} 

else 

get_state_amplitude(v) ; 

} 

/* amj #else */ 
else { /* amj */ 

_A_info[_N_atom] [4]=cos(PI/2.*_v0/v) ; /* ground state */ 
_A_info[_N_atom] [5] =-sin(PI/2 . *_v0/v) ; /* excited state */ 
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} /*amj */ 



if (0) 

{ printf ("_A_info[_N_atom] [4]=°/,g\n", _A_inf o [_N_atom] [4] ); 

printf ("_A_info[_N_atom] [5]=°/,g\n", _A_inf o [_N_atom] [5] ); 

exit(0) ; 
} 



/* CFY: change sign on sin and cos to agree with pumping-d.c */ 
/* #endif */ 

excited_prob=sqr (_A_inf o [_N_atom] [5] ) ; 
if (excited_prob>99) { 

_N_atom — ; /* throw away this new atom */ 
goto LAST.LINE; 

} /* amj add this block; for d state decayexcited_prob=100 */ 

if (GAUSSIAN.Y) y=get_height () ; 
else y=0 . ; 

_A_info[_N_atom] [0]=v; /* in u_th */ 

_A_inf o [_N_atom] [l]=t; /* arrival time in l/2g */ 

_A_inf o [_N_atom] [2] =t+INTERACTION_LENGTH*_wO/ (v*_u_th) 
*2 . *_g_over_2pi*2 . *PI ; 
/* in l/2g */ 
_A_info[_N_atom] [3]=1. ; 
if (STANDING. WAVE) 

_A_info[_N_atom] [3] *=sin((drand()-0.5)*PI) ; 
if (GAUSSIAN.Y) { 

if (sqr(y)>10.) _A_inf o [_N_atom] [3] =0 . ; 
else _A_inf o [_N_atom] [3] *=exp(-sqr (y) ) ; 
} 



_A_inf o [_N_atom] [6] =y ; /* vertical position */ 

/* Construct new coefficients _C[] */ 
max_m=two_to (_N_atom-l) *_max_photon_size ; 
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for (m=0; m<max_m; ++m) { 

_C [m+max_m] =_C [m] *_A_inf o [_N_atom] [5] ; 

_C[m]*=_A_info[_N_atom] [4] ; 

_D [m+max_m] =_D [m] *_A_inf o [_N_atom] [5] ; 

_D[m]*=_A_info[_N_atom] [4] ; 

/* 

printf ( "_C [°/„ld] =°/„lf _C [°/ ld] =°/„lf \n" , 
m, _C[m], m+max.m, _C [m+max_m] ) ; 
*/ 
} 

LAST.LINE: 

return excited_prob; 
} 



/* 

double get_arrival_time_old() 
{ 

double x, y, yy; 
do { 

x=dr and ( ) *DELTA_T_MAX ; 
y=drand() ; 

if (x<10.) yy=exp(-x); 

else yy=0. ; 

} while (y>yy) ; 

return x*_Delta_t ; 

} 

*/ 



double get_velocity() 
{ 

double x, x2, y, yy; 

if (USE_MAXWELL_BOLTZMANN) do { 
x= ( V_MAX-V_MIN) *drand ( ) +V.MIN ; 
y=drand() ; 
x2=sqr (x) ; 

if (x2<10.) yy=exp(-x2); 
else yy=0. ; 
yy*=A*(x*x2) ; 
} while (y>yy) ; 

else do { 
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x= (V_MAX-V_MIN) *drand () +V.MIN ; 
y=drand ( ) ; 
x2=sqr (x) ; 

yy=A*x*x2*exp (-x2) / ( 1+sqr (2* (x-_vO) /_delta_v) ) ; 
} while (y>yy) ; /* amj changes this block */ 

return x; 
} 



long two_to(int n) 
{ 

long x=lL; 
int i ; 

if (n<0) { 

fprintf (stderr , "Error in two_to()\n") ; 

exit(l) ; 

} 

for (i=0; i<n; ++i) x *=2L; 
return x; 
#undef V_MIN 
#undef V.MAX 



B.3 exitatom2-d.c 

void exchange_atom_with_last_one(int p) 
{ 

long max_k, mask, k, new_k, new_mO, old_mO; 
int j , bits [MAX_NUM_ATOM] ; 
double C_tmp, D_tmp; 

/* exchange p-th atom with the last atom */ 

if (p<_N_atom) { 

/* 

printf ("First exchanging °/ d-th atom with °/„d-th one\n" , 

p, _N_atom) ; 

*/ 

/* re-ordering of _C[]*/ 
max_k=two_to(_N_atom-l) ; 
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for (k=0; k<max_k; ++k) { 

/* evaluate bit patter of k */ 

mask=k; 

for (j=l; j<=_N_atom; ++j) { 
bits [ j ] =mask°/ 2 ; 

/* j-th atom state: 0=ground l=excited */ 

mask/=2 ; 

} 

/* 

printf ("k=%ld — > ", k) ; 

for (j=l; j<=_N_atom; ++j) { 

printf ("bits [°/„d]=°/„d " , j , bits[j]); 

> 

printf ("\n") ; 
*/ 



/* exchange bits[p] with bits [_N_atom] =0 */ 

bits [_N_atom] =bits [p] ; 

bits[p]=0; 

/* evaluate new k */ 
new_k=0 ; 

for (j=l; j<=_N_atom; ++j) new_k+=bits [j] *two_to(j-l) ; 
/* 

printf ("new_k=°/ ld — > ", new_k) ; 
for (j=l; j<=_N_atom; ++j) { 
printf ("bits [°/.d]=°/„d \ j, bits[j]); 
> 

printf ("\n") ; 
*/ 

/* rename indices of _C[] */ 
if ( ! (new_k==k) ) { 
new_mO=new_k*_max_photon_size ; 
old_mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) {. 

C_tmp=_C [new_m0+j] ; 

_C [new_m0+ j ] =_C [old_m0+ j ] ; 

_C[old_mO+j]=C_tmp; 

D_tmp=_D [new_m0+j] ; 

_D [new_m0+ j ] =_D [old_m0+ j ] ; 
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_D [old_mO+ j ] =D_tmp ; 
> 

/* 

printf ( "_C [°/,ld] =°/„lf _C [°/„ld] =*/.lf \n\n" , 
new_mO, _C[new_mO], old_mO, _C [old_mO] ) ; 
*/ 
} 

else { 
/* 

printf ("Same k and new_k. Not exchanged. \n\n" ) ; 

*/ 

> 

} 

/* re-order _A_inf o [] [] */ 

for (j=0; j<A_INFO_ITEMS; ++j) 

_A_info[p] [j] =_A_inf o [_N_atom] [j] ; 

} 

} 



int remove_atom(int p) 
{ 

double AO; 

long m, max_m; 

int exiting_excited=0; 

double normalize_C() ; 

exchange_atom_with_last_one(p) ; 

/* decide whether p-th atom, which is now _N_atom-th atom, 
leaves the cavity in excited or ground state */ 

/* If in ground state, do nothing. If in excited state, 
take last half of _C[] */ 
if (is_excited(&AO)) { 

max_m=two_to (_N_atom-l) *_max_photon_size ; 
for (m=0; m<max_m; ++m) { 
_C [m] =_C [m+max_m] ; 
_D [m] =_D [m+max_m] ; 

/* if ( (m°/ (_max_photon_size) )==0) */ 
/* 

printf ( "_C [°/„ld] =%lf _C ['/.Id] =°/„lf \n" , m, _C [m] , m+max_m, _C [m+max_m] ) ; 
*/ 
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} 



/*printf ("exiting in excited state (P_ex=°/ .31f ) with 
resulting _C[] array\n\n" , 
l.-AO) ;*/ 

exiting_excited=l ; 
} 

/* reduce _N_atom by 1 */ 
_N_atom — ; 

/* normalize _C[] */ 
normalize_C() ; 

/*printf ("exiting_excited=°/ d\n" , exiting_excited) ;*/ 

return exiting_excited; 

> 



int is_excited(double *A0) 
{ 

long max_m, m; 
double X; 

max_m=two_to (_N_atom-l) *_max_photon_size ; 
*A0=0 . ; 

for (m=0; m<max_m; ++m) *A0+=(sqr(_C[m] )+sqr(_D[m] )) ; 
/* ground state probability */ 
X=drand() ; 

/*printf ("A0=°/„f , X=°/„f \n" , *A0, X);*/ 
if (X> *A0) return 1; 
else { 

/*printf ("Exiting in ground state\n");*/ 

return 0; 

} 

} 



double normalize_C() 
{ 

long m, max_m; 
double norm; 
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max_m=two_to (_N_atom) *_max_photon_size ; 
norm=0 . ; 

for (m=0; m<max_m; ++m) norm+=(sqr(_C [m] )+sqr (_D [m] ) ) ; 

for (m=0; m<max_m; ++m) { 

_C [m] =_C [m] /sqrt (norm) ; 

_D [m] =_D [m] /sqrt (norm) ; 

} 

return norm; 
> 



B.4 decay2-d.c 

int cavity_decay_occurs (double dt, double *n_mean, double *n2_mean) 
{ 

double mean, sqrt_j ; 
long k, max_k; 
int j; 

double mult_f actor; 

double calculate_n_mean() ; 

double calculate_n2_mean() ; /* amj */ 

int does_cavity_decay (double , double); 

/* calculate n_mean */ 

*n_mean=calculate_n_mean() ; 

*n2_mean=calculate_n2_mean() ; 

max_k=two_to (_N_atom) ; 

if (does_cavity_decay (*n_mean, dt)) { 

/* collapse wave function */ 

for (j=l; j<_max_photon_size; ++j) { 

sqrt_j =sqrt( (double) j) ; 

for (k=0; k<max_k; ++k) { 

_C [k*_max_photon_size+ j -1] 

=_C [k*_max_photon_size+j] *sqrt_j ; 

_D [k*_max_photon_size+ j -1] 

=_D [k*_max_photon_size+j] *sqrt_j ; 

} 

} 

for (k=l; k<=max_k; ++k) { 
_C [k*_max_photon_size-l] =0 . ; 
_D [k*_max_photon_size-l] =0 . ; 
} 
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/* normalize it */ 
mean=normalize_C() ; 

/* 

printf ("Cavity decay occurred from n_mean= / f \n\n" , 

n_mean) ; 

*/ 

return 1 ; 
} 

else { 

/* decrease non-zero photon amplitudes */ 
for (j=0; j<_max_photon_size; ++j) { 

mult_f actor=l . - (double) j *_Gc_over_2pi/_g_over_2pi/4 . *dt ; 

for (k=0; k<max_k; ++k) { 

_C [k*_max_photon_size+j] *=mult_f actor ; 

_D [k*_max_photon_size+j] *=mult_f actor ; 

} 

} 

mean=normalize_C() ; 
return 0; 

> 
} 

double calculate_n_mean() 

/* calculate <n> mean photon number as well as _Pn_dist [] 

photon number distribution */ 

{ 

double mean; 
long k, max_k; 
int j; 

mean=0 . ; 

max_k=two_to (_N_atom) ; 

for (j=0; j<_max_photon_size; ++j) { 

_Pn_dist [j]=0. ; 

for (k=0; k<max_k; ++k) 

_Pn_dist [j] +=(sqr (_C [k*_max_photon_size+j] ) 
+sqr (_D [k*_max_photon_size+ j ] ) ) ; 
_Pn_dist_mean [j] +=_Pn_dist [j] ; 
mean+=_Pn_dist [j] * (double) j ; 
> 

return mean; 
> 

/* amj add the fununction calculate_n2_mean() */ 
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double calculate_n2_mean() 
/* calculate <n*n> */ 
{ 

double mean2; 
int j; 

mean2=0 . ; 

for (j=0; j<_max_photon_size; ++j) 
mean2+=_Pn_dist [j] *( (double) j )*( (double) j ) ; 

return mean2; 
> 



#define SMALL_PROB 0.3 

int does_cavity_decay (double mean, double dt) 
{ 

double prob; 

/* calculate decay rate */ 

prob=_Gc_over_2pi/_g_over_2pi/2 . *dt*mean ; 

if ((prob > 0.3)&&(prob < 10.)) prob=l . -exp(-prob) ; 

if (prob>10.) prob=l. ; 

if (drand()<prob) return 1; 

else return 0; 

> 



B.5 evolve2-2-d.c 

void evolve_wave_f unction(double dt, double t) 

/* Runge-Kutta 4th order algorithm */ 

{ 

double dC, dD, x, y, g [MAX_NUM_AT0M] ; 
long mO, mOj , k_new, k, mask; 
int j , i ; 

double xp, tO, omega [MAX_NUM_ ATOM] ; 

for (k=0; k<_max_k; ++k) { 

/* evaluate bit pattern of k */ 

mask=k; 
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_detuning_sum [k] =0 . ; /* freq/2g unit */ 
for (j=l; j<=_N_atom; ++j) { 
_bits[k] [ j ] =mask°/ 2 ; 

/* j-th atom state: 0=ground l=excited */ 
if (_bits[k] [j]>0) 
_detuning_sum [k] += (_CA_detuning- 
_A_info[j] [0] *_mean_shif t_over_2g/_v0) ; 
else 

_detuning_sum [k] -= (_CA_detuning- 
_A_info[j] [0] *_mean_shif t_over_2g/_v0) ; 
mask/=2; 

} /* amj add _v0 jun 16,98 */ 

_detuning_sum [k] *=0 . 5 ; 

} 

/* calculate coupling constant of i-th atom */ 
for (i=l; i<=_N_atom; ++i) { 

g[i] =_A_inf o [i] [3] ; /* of i-th atom, already including 
exp(-sqr(y)) factor */ 

x=(2 . *(t-_A_inf o [i] [1] ) /* calculate its x position at t 
/(_A_info[i] [2]-_A_info[i] [1])-1.) 
*INTERACTI0N_LENGTH/2. ; 

/* -INTERACT I 0N_LENGTH/ 2 < x < + INTERACT I ON.LENGTH/ 2 */ 
if (GAUSSIAN.X) { 

if (sqr(x)<10.) g [i] *=exp(-sqr (x) ) ; 

else g[i]=0. ; 

} 

else { /* Top hat in X */ 

if (f abs (x)>INTERACTI0N_LENGTH/2 . ) g[i]=0.; 
> 

/* 

if(i==l ) printf ("t=°/„g: g[l]=°/„g\n", t, g[l]); 
*/ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 
/* check if the atom is still in the pump beam */ 
if (fabs(x+PUMP_P0S)<PUMPING_LENGTH/2.) { 
xp= (x+PUMP_P0S) /PUMP.WAIST ; 

if (sqr(xp)<10. ) omega [i] =_omega_Rabi*exp (-sqr (xp) ) ; 

else omega [i] =0 . ; 
} /* atom is still in the pump */ 
else omega [i] =0 . ; 
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} /* the pump overlaps with the cavity */ 
} 

/* First, calculate yl=yO+y ' (t ,y0)*dt/2 for RK4 */ 
/* or just y=yO+y J (t,yO)*dt for non-RK4 */ 

for (k=0; k<_max_k; ++k) { 
mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) { 

/* evaluate two atom-field interaction terms */ 
/* first, initialize with detuning term */ 
mOj=mO+j ; 

dC=-_D [mOj] *_detuning_sum[k] ; /* detuning in 2g unit */ 
dD=_C [mO j ] *_detuning_sum [k] ; 

if (j > 0) { /* a S- term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==0) { 

/* only when i-th atom in ground state */ 
k_new=k+two_to (i-1) ; 

/* H=hbar*g*(a+S- + aS+) used here */ 

dC+=_D [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j ) *0 . 5 ; 
dD-=_C [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D [k_new*_max_photon_size+j] *0 . 5; 

dD-=omega[i] *_C [k_new*_max_photon_size+j] *0 . 5; 

} /* the pump overlaps with the cavity */ 

} /* if i-th atom in ground state */ 

} /* for i loop */ 

} /* if */ 

if (j<_max_photon_size-l) { /* a S+ term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==D { 

/* only when i-th atom in excited state */ 
k_new=k-two_to (i-1) ; 

dC+=_D [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) ( j +1) ) *0 . 5 ; 
dD-=_C [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) ( j +1) ) *0 . 5 ; 
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/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D [k_new*_max_photon_size+j] *0 . 5; 

dD-=omega[i] *_C [k_new*_max_photon_size+j] *0 . 5; 

} /* the pump overlaps with the cavity */ 

} /* if i-th atom in excited state */ 

} /* for i loop */ 

} /* if */ 

dC*=dt; /* dt in l/2g unit */ 
dD*=dt; /* dt in l/2g unit */ 
#ifdef USE_RK4_ALG0RITHM 
dC/=2 . ; 
dD/=2 . ; 
#endif 

_C_1 [m0j]=_C[m0j]+dC; /* dt in l/2g unit */ 
_D_1 [m0j]=_D[m0j]+dD; /* dt in l/2g unit */ 
} /* for j loop */ 
} /* for k loop */ 



#ifndef USE_RK4_ALG0RITHM /* not USE_RK4_ALG0RITHM */ 

for (k=0; k<_max_k; ++k) { 

mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) { 

mOj=mO+j ; 

_C[m0j]=_C_l[m0j] ; 
_D[m0j]=_D_l[m0j] ; 
> 
} 

#else /* do USE_RK4_ALG0RITHM */ 

/* Second, calculate y2=y0+y' (t+dt/2,yl)*dt/2 */ 
t+=dt/2; 

for (k=0; k<_max_k; ++k) { 
mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) { 
mOj=mO+j ; 

/* evaluate two atom-field interaction terms */ 

/* first, initialize with detuning term */ 

dC=-_D_l [mOj] *_detuning_sum[k] ; /* detuning in 2g unit */ 

dD=_C_l [mOj] *_detuning_sum[k] ; 
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if (j > 0) { /* a S- term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==0) { 

/* only when i-th atom in ground state */ 
k_new=k+two_to (i-1) ; 

dC+=_D_l [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j )*0 .5; 
dD-=_C_l [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D_1 [k_new*_max_photon_size+j] *0 . 5 ; 

dD-=omega[i] *_C_1 [k_new*_max_photon_size+j] *0 . 5 ; 

} /* the pump overlaps with the cavity */ 

} /* if i-th atom in ground state */ 

}/* for loop */ 

} 

if (j<_max_photon_size-l) { /* a S+ term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==D { 

/* only when i-th atom in excited state */ 
k_new=k-two_to (i-1) ; 

dC+=_D_l [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 
dD-=_C_l [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D_1 [k_new*_max_photon_size+j] *0 . 5 ; 

dD-=omega[i] *_C_1 [k_new*_max_photon_size+j] *0 . 5 ; 

} /* atom is still in the pump */ 

} /* if i-th atom in excited state */ 

} 

} 

_C_2[m0j]=_C[m0j]+dC*dt/2. ; /* dt in l/2g unit */ 
_D_2[m0j]=_D[m0j]+dD*dt/2. ; /* dt in l/2g unit */ 
} /* for j loop */ 
}/* for k loop */ 

/* Third, calculate y3=y0+y' (t+dt/2,y2)*dt 



212 



for (k=0; k<_max_k; ++k) { 
mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) { 
mOj=mO+j ; 

/* evaluate two atom-field interaction terms */ 

/* first, initialize with detuning term */ 

dC=-_D_2 [mOj] *_detuning_sum[k] ; /* detuning in 2g unit */ 

dD=_C_2[m0j] *_detuning_sum[k] ; 

if (j > 0) { /* a S- term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==0) { 

/* only when i-th atom in ground state */ 
k_new=k+two_to (i-1) ; 

dC+=_D_2 [k_new*_max_photon_size+j -1] *g [i] *sqrt ( (double) j ) *0 . 5 ; 
dD-=_C_2 [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D_2 [k_new*_max_photon_size+j] *0 . 5 ; 

dD-=omega[i] *_C_2 [k_new*_max_photon_size+j] *0 . 5 ; 

} /* atom is still in the pump */ 

} /* if i-th atom in ground state */ 

} /* for loop */ 

> 

if (j<_max_photon_size-l) { /* a S+ term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==D { 

/* only when i-th atom in excited state */ 
k_new=k-two_to (i-1) ; 

dC+=_D_2 [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 
dD-=_C_2 [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D_2 [k_new*_max_photon_size+j] *0 . 5 ; 

dD-=omega[i] *_C_2 [k_new*_max_photon_size+j] *0 . 5 ; 

} /* atom is still in the pump */ 

} /* if i-th atom in excited state */ 

> 

} 

_C_3[m0j]=_C[m0j]+dC*dt; /* dt in l/2g unit */ 
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_D_3[m0j]=_D[m0j]+dD*dt; /* dt in l/2g unit */ 
} /* for j loop */ 
} /* for k loop */ 

/* Fourth, calculate yf =(yl+2*y2+y3-y0)/3+y' (t+dt,y3)*dt/6 */ 
t+=dt/2; 

for (k=0; k<_max_k; ++k) { 
mO=k*_max_photon_size ; 

for (j=0; j<_max_photon_size; ++j) { 
mOj=mO+j ; 

/* evaluate two atom-field interaction terms */ 

/* first, initialize with detuning term */ 

dC=-_D_3 [mOj] *_detuning_sum[k] ; /* detuning in 2g unit */ 

dD=_C_3 [mOj] *_detuning_sum[k] ; 

if (j > 0) { /* a S- term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==0) { 

/* only when i-th atom in ground state */ 
k_new=k+two_to (i-1) ; 

dC+=_D_3 [k_new*_max_photon_size+j-l] *g [i] *sqrt ( (double) j )*0 .5; 
dD-=_C_3 [k_new*_max_photon_size+j -1] *g [i] *sqrt ( (double) j ) *0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 

dC+=omega[i] *_D_3 [k_new*_max_photon_size+j] *0 . 5 ; 

dD-=omega[i] *_C_3 [k_new*_max_photon_size+j] *0 . 5 ; 

} /* atom is still in the pump */ 

} /* if i-th atom in ground state */ 

} /* for loop */ 

} 

if (j<_max_photon_size-l) { /* a S+ term */ 
for (i=l; i<=_N_atom; ++i) { 
if (_bits[k] [i]==D { 

/* only when i-th atom in excited state */ 
k_new=k-two_to (i-1) ; 

dC+=_D_3 [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 
dD-=_C_3 [k_new*_max_photon_size+j+l] *g [i] *sqrt ( (double) (j+1) )*0 . 5 ; 

/* include pump beam effect on lasing */ 

if ( ( PUMP I NG_ S I MUL ATED ) && (PUMPING_EFFECT_ON) ) { 
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dC+=omega[i] *_D_3 [k_new*_max_photon_size+j] *0 . 5 ; 
dD-=omega[i] *_C_3 [k_new*_max_photon_size+j] *0 . 5 ; 
} /* atom is still in the pump */ 
} /* if i-th atom in excited state */ 
} 
} 

_C [mOj ] = (_C_ 1 [mO j ] +_C_2 [mO j ] *2 . +_C_3 [mO j ] 
/* dt in l/2g unit */ 

_D [mO j ] = (_D_ 1 [mO j ] +_D_2 [mO j ] *2 . +_D_3 [mO j ] 
/* dt in l/2g unit */ 
} /* for j loop */ 
} /* for k loop */ 

t-=dt; /* unchanged t */ 

#endif 
} 



B.6 spont2-d.c 

int decide_atom_decay (int *does_it_decay , double dt) 
{ 

int i , n=0 ; 
double prob; 

for (i=l; i<=_N_atom; ++i) { 

prob=_Ga_over_2g*dt*_population [i] ; 

if ((prob>0.3)&&(prob<10. )) prob=l . -exp(-prob) ; 

if (prob>10.) prob=l.; 

if (drand()<prob) { 

does_it_decay [i] =1 ; 

n++; 

/* 

printf ("%d-th atom undergoes spontaneous emission. \n" , i) ; 

*/ 

} 

else does_it_decay[i]=0; 
} 

return n; 
} 



-_C[m0j])/3.+dC*dt/6. ; 
-_D [m0j])/3.+dD*dt/6. ; 
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int how_many_atoms_decay (double dt) 
{ 

long max_k, k, mask, new_k, m, m2; 

int bits [MAX_NUM_ATOM] , j, i, does_it_decay [MAX_NUM_ATOM] ; 
double mean; 
int num_photons; 

/* initialize population array */ 

for (i=0; i<MAX_NUM_ATOM ; ++i) .population [i] =0 . ; 

max_k=two_to(_N_atom) ; 

/* calculate excited state population */ 

for (k=0; k<max_k; ++k) { 

mask=k; 

for (i=l; i<=_N_atom; ++i) { 
/* evaluate i-th bit pattern */ 
bits [i]=mask°/ 2; 
mask/=2; 

if (bits[i]==l) { 

for (j=0; j<_max_photon_size; ++j) { 

m=k*_max_photon_size+j ; 

.population [i] +=(sqr(_C [m] )+sqr(_D [m] ) ) ; 

} 

} 

} 

} 

num_photons=decide_atom_decay (does_it_decay , dt) ; 

for (k=0; k<max_k; ++k) { 
mask=k; 

for (i=l; i<=_N_atom; ++i) { 
/* evaluate i-th bit pattern */ 
bits [i]=mask°/ 2; 
mask/=2 ; 

if (does_it_decay [i] ==1) { 
/* collapse wave function */ 
if (bits[i]==l) { 
new_k=k-two_to(i-l) ; /* exist */ 
for (j=0; j<_max_photon_size; ++j) {. 
m2=new_k*_max_photon_size+j ; 
m=k*_max_photon_size+j ; 
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_C [m2] +=_C [m] ; /* down shift */ 
_C [m] =0 . ; 

_D [m2] +=_D [m] ; /* down shift */ 
_D [m] =0 . ; 

} /* otherwise no new term */ 
} 

else 

for (j=0; j<_max_photon_size; ++j) { 

_C [k*_max_photon_size+ j ] =0 . ; 

_D [k*_max_photon_size+ j ] =0 . ; 

> 

} 

else {/* does not decay */ 
if (bits[i]==l) { 

/* decrease excited state amplitude */ 
for (j=0; j<_max_photon_size; ++j) { 
m=k*_max_photon_size+j ; 
/* 

printf ( M C[ /„ld]=°/„f — >", m, _C[m]); 
*/ 

_C[m]*=(l.-_Ga_over_2g/2.*dt) ; 
_D[m]*=(l.-_Ga_over_2g/2.*dt) ; 
/* 

printf ("C[y„ld]=°/„f\n", m, _C[m]); 

*/ 

} 

} 

else { 

/* do nothing */ 
/* 

for (j=0; j<_max_photon_size; ++j) { 

m=k*_max_photon_size+j ; 

printf ("C[ /old]=°/.f\n" ) m, _C[m]); 

> 

*/ 

} 

} 

} 

> 

mean=normalize_C() ; 
return num_photons ; 
> 
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B.7 pumping-d.c 



/* This part incorporates coherent pumping process */ 

/* #define SELF_TESTING this part removed by CFY because it 

/* did not compile correctly - old version still exists in pumping-d.old.c */ 

int TEST_PUMP=0 ; 

double FUDGE_FACT0R=1 ; 

void get_state_amplitude (double v) 
{ 

double dt, Dt, a, b, t, tl, t2, tp, xl, x2, 
Ga_over_Omega_R, prob; 

int i ; 

char *output_name = "pump. out"; 
FILE *fp; 

/* integration from -PUMPING.LENGTH/2 to PUMP_P0S-INTERACTI0N_LENGTH/2 */ 

Ga_over_0mega_R=_Ga_over_2g*_t_int*PUMP_WAIST/PI; 

xl=-PUMPING_LENGTH/2. ; 

x2=PUMP_P0S-INTERACTI0N_LENGTH/2. ; /* amj /2 */ 

/* xl, x2 in wO unit */ 

/* printfC'*** v=°/„g, _vO=°/ g\n", v, _v0) ; */ 



FUDGE.FACTOR = 1.187; 

tp=sqrt(PI)/(v/_vO) * FUDGE.FACTOR; 

tl=xl/PUMP_WAIST*tp; 

t2=x2/PUMP_WAIST*tp ; 

/* tp, tl, t2 in l/0mega_R unit */ 

t=tl; 
dt=0.01; 

a=l . ; /* ground state amplitude */ 
b=0.; /* excited state amplitude */ 

do { 

Dt=dt*exp(-sqr(t/tp)) ; 
a+=b/2 . *Dt ; 
b+=-a/2 . *Dt ; 
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t+=dt ; 



/* amj test starts jul 27, 98 */ 

if (TEST.PUMP) 
{ 

for (i=0; i<2; ++i) 
{ 

if (i==0) 

fp=f open ( output _name, "a"); 
else fp=stderr; 

fprintf(fp, "°/„lf °/„lf °/„lf \n", t, a, b) ; 
if(i==0) fclose(fp); 

} 
} 



/* amj test ends jul 27, 98 */ 



/* calculate decay probability */ 
/* amj #ifdef INCLUDE_ATOM_DECAY */ 
if ( INCLUDE_ATOM_DECAY) { 

prob=Ga_over_Omega_R*dt*sqr (b) ; 

if ((prob>0.3)&&(prob<10. )) prob=l . -exp(-prob) ; 

if (prob>10.) prob=l.; 

/* generate a random number */ 

if (drand()<prob) { 

/* atom decays */ 

b=0. ; 

a=l.; 

}} 

if ( INCLUDE_D_STATE_DECAY) { 

prob=D_STATE_DECAY_CORRECTION*Ga_over_Omega_R*dt*sqr(b) ; 
if ((prob>0.3)&&(prob<10. )) prob=l . -exp(-prob) ; 
if (prob>10.) prob=l. ; 

/* generate a random number */ 

if (drand()<prob) b=10; /* amj by definition b=10 for d state decay */ 
} /* amj adds this block */ 
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/* amj else O */ /* if not decaying, do nothing */ 
/* amj #endif */ 

} while (t<t2) ; 

/* amj test starts jul 27, 98 */ 
/* 

for (i=0; i<2; ++i) 
{ 

if (i==0) 

fp=f open ("pumping. out" , "a") ; 
else fp=stderr; 

fprintf(fp, "after this point atom in interaction region\n"); 
if (i==0) fclose(fp); 

} 

*/ 

if (TEST.PUMP) 
exit(O) ; 

/* amj test ends jul 27, 98 */ 

_A_info[_N_atom] [4]=a; 
_A_inf o [_N_atom] [5]=b; 



} 



B.8 d_state2-d.c 

int decide_d_state_decay(int *does_it_decay , double dt) 
{ 

int i , n=0 ; 
double prob; 

for (i=l; i<=_N_atom; ++i) { 

/* calculate decay probability to D state */ 

prob=_Ga_over_2g*dt*_population [i] *D_STATE_DECAY_CORRECTION ; 

if ((prob>0.3)&&(prob<10. )) prob=l . -exp(-prob) ; 

if (prob>10.) prob=l.; 

if (drand()<prob) { 

does_it_decay [i] =1 ; 
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n++; 
/* 

printf ("°/ d-th atom undergoes d-state decay. \n" , i) ; 

*/ 

> 

else does_it_decay[i]=0; 
} 

return n; 
} 

int how_many_d_state_decay (double dt) 
{ 

long max_k, k, mask, m, max_m; 

int bits.i, j, i, does_it_decay [MAX_NUM_ATOM] ; 

double mean; 

int num_decays; 

/* initialize population array */ 

for (i=0; i<MAX_NUM_ATOM ; ++i) .population [i] =0 . ; 

max_k=two_to (_N_atom) ; 

/* calculate excited state population */ 

for (k=0; k<max_k; ++k) { 

mask=k; 

for (i=l; i<=_N_atom; ++i) { 
/* evaluate i-th bit pattern */ 
bits_i=mask°/ 2 ; 
mask/=2; 

if (bits_i==l) { 

for (j=0; j<_max_photon_size; ++j) { 

m=k*_max_photon_size+j ; 

.population [i] +=(sqr(_C [m] )+sqr(_D [m] ) ) ; 

} 

} 

} 

} 

num_decays=decide_d_state_decay (does_it_decay , dt) ; 

for (i=l; i<=_N_atom; ++i) { 
/* evaluate i-th bit pattern */ 
if (does_it_decay [i] ==1) { 

/* i-th atom decays to d-state — > i-th atom should be in the 
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excited state just before the decay. When D-state decay occurs 
for i-th atom, first we exchange the atom with the last one. 
Now the i-th atom is labeled as the last one. Since the new last 
one is in D-state, it is as if it were not in the cavity. We 
should remove the atom from the wave function. Just before the 
removal, the atom should be in the excited state. This is equivalent 
to the atom exciting cavity in the excited state prematurally . */ 

exchange_atom_with_last_one(i) ; 

/* the last atom must be in the excited state */ 

max_m=two_to (_N_atom-l) *_max_photon_size ; 

for (m=0; m<max_m; ++m) { 

_C[m]=_C[m+max_m] ; 

_D [m] =_D [m+max_m] ; 

> 

_N_atom — ; 
> 

else {/* does not decay */ 
max_k=two_to (_N_atom-l) ; 
for (k=0; k<max_k; ++k) { 
mask=k; 

/* find the i-th bit pattern of k */ 
for (j=l; j<=i; ++j) { 
bits_i=mask°/ 2 ; 
mask/=2 ; 
} 

if (bits_i==l) { 

for (j=0; j<_max_photon_size; ++j) { 

m=k*_max_photon_size+j ; 

/* decrease excited state amplitude */ 

_C [m] *=(1 . -_Ga_over_2g/2 . *dt*D_STATE_DECAY_CORRECTION) ; 

_D [m] *=(1 . -_Ga_over_2g/2 . *dt*D_STATE_DECAY_CORRECTION) ; 

} /* for j loop */ 

} /* if */ 

} /* for k loop */ 

} /* if */ 

} /* for i loop */ 

mean=normalize_C() ; 
return num_decays; 
> 
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B.9 Sample input file 



50.000000 # Ga_over_2pi 
791.000000 # lambda 
10.000000 # rO 
10 # finesse 
0.100000 # L 

200 3 # _max_photon_size _max_atoms 

0.05000 # dtO 

10000 # tot_ atom. count 

# DEBUG.LEVEL 

INCLUDE_ATOM_DECAY 

# INCLUDE_D_STATE_DECAY 

# STANDING. WAVE 
10.000000 # tilt_angle 

1 # num_detuning 
1.000000 # detuning [0] 
4.000000 # INTERACTION_LENGTH 

# CIRCULAR.BEAM 
0.0001 # REC_HALF_HEIGHT 

0.443 # effective_V_to_apparent_V 

1 # num.N.eff 

.1 # N_ef f [0] *** Don't forget summary filename below *** 

1 # num_u_th 

400.000000 # u_th[0] 

# M0N0.VEL0CITY 

# USE_MAXWELL_B0LTZMANN 

0.100000 # _delta_v 

1.000000 # _v0 

# GENERATE_PH0T0N_FILE 

# PUMP I NG_ S I MUL ATED 

run40.data # summary filename 

Note: this program may not treat velocity distributions and detuning correctly; 
in this thesis simulation results were only cited for the monovelocity, resonant case. 
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Appendix C 



VPascal sequencing program 



C.l Cavity scanning experiment 



{ microlaser_experiment . v } 
{ C. Fang-Yen } 

{ Variables } 
var 

Zl ; plotimage ; bgimage ; img; 

PlotO ; Plotl ; Plot2 ; Plot3 ; expo sure .time ; 

imagename; backgroundname ; scatteringname ; scatteringimage ; 

rowplotname ; colplotname ; 

do_dataplot ; dataplotname ; crop_colplot ; 

i; index; x; { scanned voltage, etc. } 

x_start; x_stop; x_increment; x_use; 

micp; mirp; mox; moy; 

writedata; textwin; textwinname; 

do_serial; total_counts_x; total_counts_y ; num_pixels; 
wherexmax; whereymax; iylimit; ixlimit; sumO; suml; sum2; 
rezero; rezero_sample ; moymox; y_third_max_minus ; y_third_max_plus ; 
ii; do_gaussian_f it ; closest; fitmin; fitmax; 
sigma; logmox; logmoxregion; a; line; linelOO; 
test; absmox; data_x; data_y; data_x_temp ; data_y_temp ; 
doing_background_now; signal; previous_signal ; ratio; serial_error ; 
serial_error_msg; dataline; current_position; background_plot ; 
n_periods; data_array; n_frames; data_max; data_bkgnd; where_max; 
data_bkgnd_sample ; v_t_553_on; v_t_791_on; v_t_553_off; v_t_791_off; 
use_background; plotimage_no553 ; dataf ilename ; do_blackbody ; do_microlaser ; 
save_images ; image_array; image_array_f ilename ; cts_per_ms; 
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procedure take_picture(roi_index; var img) ; 
begin 

{ Set Gain Index 3 } 
pvcSetGain( 3 ) ; 
exposure_time := 300; 
pvcSetExpTime ( exposure_time) ; 
roi_index := 3; { override } 

if(roi_index = 0) then 

img := pvcSequence( 1,386,228,423,257,1,1 ) ; {small} 
if(roi_index = 1) then 

img := pvcSequence( 1,304,169,501,323,1,1 ) ; {medium} 
if(roi_index = 2) then 

img := pvcSequence( 1,0,0,767,511,1,1 ); { full CCD } 

if(roi_index = 3) then 

img := pvcSequence( 1,352,236,454,270,1,1 ); { small for 
slit} 

end; 
begin 

-[ ******************* INITIALIZATION **************** } 

datafilename := 'data.tif; { 'data.tif = generic } 
image_array_f ilename := ' image .tif ' ; 

do_blackbody := False; { Record location of blackbody peak in col plot } 
use_background := False; { else use scattering subtraction } 
do_microlaser := True; { else skip cavity scan } 
save_images := False; { save fluorescence images } 

imagename : = ' image ' ; 
rowplotname := 'rowplot'; 
colplotname := 'colplot'; 
dataplotname := 'dataplot'; 
crop_colplot := 0; 
backgroundname := 'background'; 
scatteringname := 'scattering'; 
index := 0; 
x_use := True; 
x.start :=1450; 
x_stop := 1500; 
x_increment := 2; 
x := x_start; 
writedata := True; 
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textwinname : = ' output ' ; 
do_serial := True; 
rezero := False; 
rezero_sample := 10; 
do_gaussian_f it := False; 
data.x := MakeLinear (0 . , . 0, 1000) ; 
data.y := MakeLinear (0 . , . 0, 1000) ; 
do_dataplot := False; 
doing_background_now := 0; 

v_t_553_on := 0.65; { V_t 553nm on voltage } 
v_t_791_on := 1.0; { V_t 791nm on voltage } 
v_t_553_off := -0.1; { V_t 553nm off voltage } 
v_t_791_off := -0.1; { V_t 791nm off voltage } 



if (writedata) then 
begin 

if EditorExists (textwinname) then 

GetEditor (textwinname , textwin) 
else textwin := CreateEditor ( 'output'); 
WriteLn(textwin, ' ') ; 
end; 

{ ******* INITIALIZE SERIAL PORTS ****** } 

SetTxEnd( chr(13) ) ; 
SetRxEnd( chr(13) ) ; 
SetRxTimeout(5000) ; 

WriteLnC ') ; 

OpenSeriaK 1 , 19200 ,8,NoParity , 1 ) ; { SR400 photon counter A } 
serial_error := SerialError; 

if serial_error <> then 
Halt( 'Failed to open serial connection 1: error ', serial_error) 
else WriteLn( 'Serial connection 1 ok'); 

SelectPort(l) ; 
TxFlush ; 

Transmit ( chr(13)) ; 

Transmit ( chr(13)) ; 
TxFlush ; 
RxFlush; 
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Transmit (' pi 1'); 

dataline := RxString; 
WriteLn('SR400A: P0RT1=' , dataline); 

Transmit ( 'np' ) ; 

dataline := RxString; 
n_periods := val (dataline) ; 
WriteLn('SR400A: NPERIODS=', n.periods) ; 

{ 

CloseSerial(l) ; 

Halt ( " ) ; 

} 

if x_increment <> then 

n_frames := Floor ( (x_stop - x_start)/x_increment)+l 
else 

n_frames := 100; 

data_array := CreateArray (single, n_periods+l, n_frames, 1); 

OpenSeriaK 2,9600,8,NoParity, 1 ) ; { IMS stepper driver } 
serial_error := SerialError; 

if serial_error <> then 
Halt( 'Failed to open serial connection 2: error ', serial_error ) 
else WriteLn( 'Serial connection 2 ok'); 

{ ****** INITIALIZE BLOCKER ****** } 

SelectPort(2) ; 

Transmit ( ' ' ) ; 

Delay(lOO) ; 
RxFlush; 

WriteLn( 'RxWaiting=' ,RxWaiting) ; 

Transmit ( ' Z ' ) ; 

Delay(lOO) ; 

dataline := RxString; 

WriteLn( 'Current position ', dataline); 
dataline := StrParse (dataline , 'Z'); 
current_position := Val (StrStrip (dataline) ) ; 
WriteLn( 'Current position Z=', current_position) ; 
Transmit (' X' ) ; 
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dataline := RxString; 
WriteLn(dataline) ; 
dataline := RxString; 
WriteLn (dataline) ; 

{ 

CloseSerial(l) ; CloseSerial (2) ; 

Halt ( 'done' ) ; 

} 

{ ************* TAKE 1 FRAME, SET UP IMAGE AND PLOT WINDOWS ************ } 
{ Acquire data } 

{ plotimage := pvcSequence( 1,386,228,423,257,1,1 ) ; } 
take_picture(0, plotimage); 
Show( plotimage) ; 

num_pixels := GetXSize (plotimage) * GetYSize (plotimage) ; 

if save_images then image_array := CreateArray (GetXSize (plotimage) , 
GetYSize (plotimage) , n_f rames) ; 

WriteLn (FileExpand( 'background.tif ' ) ) ; 
{ Halt; } 

if use_background then 
if not ImageExists(backgroundname) then 
begin 

if Query ('No background: use this? ' )=id_Yes then 
begin 

Save (plotimage , 'background.tif ' ) ; 
OpenCbackground.tif , background_plot) ; 
Show (background_plot , ' background ' ) ; 
Delete (plotimage) ; 
end; 
Halt; 
end; 

{ 

if use_background = False then 

if not ImageExists(scatteringname) then Halt ('No scattering file'); 

} 

{ Get the active frame number } 

Zl := GetActiveFrame ( plotimage ) ; 
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{ Create a plot window } 
if PlotExists (rowplotname) then 
begin 

GetPlot (rowplotname , PlotO) ; 
{ Halt (' ok'); } 
end 

else 
begin 

PlotO := CreatePlot ( rowplotname ); 
SetTitle( PlotO, 'Row Plot' ) ; 
SetSubTitle( PlotO ,GetName ( plotimage ) ) ; 
SetXLabeK PlotO,' X Position' ) ; 
SetYLabeK PlotO, 'Mean Intensity' ) ; 
SetGrid( PlotO, gs_Both ) ; 
end; 

{ Row plot } 

SetPlotColor( col.Black ) ; 

Plot( PlotO, MeanOverY( plotimage [..,.., Zl] ) ) ; 

{ Create data plot window } 
if PlotExists (dataplotname) then 
begin 

GetPlot (dataplotname , Plot2) ; 

end 
else 

Plot2 := CreatePlot (dataplotname) ; 

{ Create a plot window } 
if PlotExists (colplotname) then 
begin 

GetPlot (colplotname, Plotl) ; 

{ Halt ('ok'); } 
end 
else begin 

Plotl := CreatePlot ( colplotname ) ; 

SetTitle( Plotl , 'Column Plot' ) ; 

SetSubTitle( Plotl,GetName( plotimage ) ) ; 

SetXLabeK Plotl, 'Y Position' ) ; 

SetYLabeK Plotl, 'Mean Intensity' ) ; 

SetGrid( Plotl ,gs_Both ) ; 
end; 
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{ Col plot } 

SetPlotColor( col_Black ) ; 

Plot( Plotl,MeanOverX( plotimage [..,.., Zl] ) ) ; 
{ *** INITIALIZE SR400 PARAMETERS *** } 
{ to be written } 



{ ******************* MAIN LOOP BEGINS ************************ } 
repeat 

SelectPort(2) ; 

Transmit (StrCatCR ', Str(x))); { SET BLOCKER POSITION } 

Delay (1000. 0*Abs (current.position - x)/400.0); { Wait for it to move } 

current_position := x; 

{ Turn on 553nm probe, turn off 791nm pump} 
SelectPort(l) ; 

Transmit(StrCat('pl 1, Str (v_t_553_on) ) ) ; 
Transmit(StrCat('pl 2, Str (v_t_791_of f ) ) ) ; 

Delay(lOOO) ; 

{ Acquire data } 
take_picture(0, plotimage); 
Show( plotimage ) ; 

{ Turn off 553nm probe, turn on 791 pump } 
SelectPort(l) ; 

Transmit (StrCatC pi 1, Str (v_t_553_of f ) ) ) ; 
Transmit (StrCatC pi 2, Str (v_t_791_on) ) ) ; 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 

{ plotimage [.., .. crop_colplot ,.. ] := mox[0,crop_colplot,0] ; } 

Get Image ( backgroundname , bgimage ) ; 
Getlmage( scatteringname , scatteringimage ); 
{ Perform an image subtraction } 
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if use_background then plotimage := plotimage - bgimage ; 

if not use_background then 
begin 

take_picture(0, plotimage_no553) ; 

plotimage := plotimage - plotimage_no553 - 0; 

end; 

{ Rezero: subtract avg of first 10 pixels in col plot } 
if rezero then begin 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 
moymox : = MeanOverY (mox [ . . , . . 10 , . . ] ) ; 
plotimage := plotimage - moymox [0,0,0] ; 
{ mox[0,crop_colplot,0] ; } 
end; 

if save_images then 
begin 

image_array [..,.., index] := plotimage; { store plotimage in array } 
if Save (image_array , image_array_f ilename) <> file_0k 
then Halt ('Save error'); 

end; 

{ Row plot } 

ClearPlot(PlotO) ; 

SetPlotColor( col_Black ) ; 

moy := MeanOverY ( plotimage [..,.. ,Z1] ); 

Plot( Plot0,moy ) ; 
mirp := MaxOf (moy) ; 

wherexmax := 0; 

ixlimit := GetXSize (moy) -1 ; 

for i := to ixlimit do 

begin 

if moy [i, 0,0] = mirp 

then wherexmax := i; 
{ WriteLn(str(i)) ; } 
end; 

sumO := 0; suml := 0; sum2 := 0; 

for i :=0 to ixlimit do 

begin 
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sumO := sumO + moy[i,0,0] ; 
suml := suml + moy[i,0,0] * i; 
sum2 := sum2 + moy[i,0,0] * i * i; 
end; 

WriteLn('x: avg=', suml/sumO, ' sigma^, 
Sqrt (sum2/sum0 - suml*suml/ (sumO*sumO) ) , ' sum=', sumO) ; 

WriteLn(x, ' ', sumO) ; 

{ WriteLn(textwin, 'max at x=', wherexmax) ; } 

{ WriteLn(textwin, MaxOf (moy [0 . . 10] ) ) ; } 
total_counts_y := MeanOf (moy) * num_pixels; 
SetTitleC PlotO, StrCatCStrCatORow Plot, max=', Str(mirp)), 
StrCatC at ', Str (wherexmax) ) ) ) ; 

{ Col plot } 

ClearPlot(Plotl) ; 

SetPlotColor( col_Black ) ; 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 

{ WriteLn(textwin, GetYSize(mox)) ; } 

mox [ . . ,0 . . crop_colplot , . . ] : = mox [ . . , crop_colplot , . . ] ; 
{ if crop_colplot>0 then 

for i := 1 to crop_colplot-l do 
begin 
WriteLn(i) ; 

mox[i] := mox [crop_colplot] ; 
end; } 

Plot( Plot 1, mox ) ; 
micp := MaxOf (mox) ; 

total_counts_x := MeanOf (mox) * num_pixels; 

whereymax := 0; 

iylimit := GetYSize (mox) -1 ; 

for i := to iylimit do 

begin 

if mox[0,i,0] = micp 

then whereymax := i; 
{ WriteLn(str(i)) ; } 
end; 

{ Find 1/3 max points } 
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closest := 10000; 
for ii := whereymax to iylimit do 
begin 

if Abs(mox[0,ii,0] - micp/3.0) < closest then 
begin 

closest := Abs(mox[0,ii,0] - micp/3.0); 
y_third_max_plus := ii; 
end; 
end; 

closest := 10000; 

for ii:=0 to whereymax do 



WriteLn('y_half _max_plus str (y_third_max_plus) ) ; 
WriteLn( ; y_half _max_minus ', str (y_third_max_minus) ) ; 



{ Gaussian fitting } 
sigma :=0; 

if do_gaussian_f it then begin 

absmox := sqrt(sqr(mox))+0.01; 
logmox := Ln (absmox) ; 
fitmin := y_third_max_minus ; 
fitmax := y_third_max_plus ; 

line := MakeLinear (fitmin -0.0 , fitmax - 0.0, fitmax - fitmin + 1) ; 
linelOO := MakeLinear (0, iylimit+0 . 0, (iylimit+1 . 0) ) ; 
a := PolyFit (line , logmox [0, fitmin. .fitmax, 0] ,2) ; 

{ WriteLn(a) ; } 
{ 

WriteLn(str (line) ) ; 

WriteLn (str (logmox [0, fitmin. . fitmax, 0] )) ; 
WriteLn(str(a[0,0,0] )) ; 
WriteLn(str(a[l,0,0] )) ; 
WriteLn (str (a [2, 0,0])) ; 
WriteLnC '); } 



begin 

if Abs(mox[0,ii,0] - micp/3 
begin 

closest := Abs(mox[0,ii,0] 
y_third_max_minus := ii; 
end; 

end; 

{ 



0) < closest then 



- micp/3.0) ; 



} 
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{ logmoxregion := logmox[0,f itmin. . f itmax,0] ; 
WriteLn(logmoxregion) ; 
WriteLn(str(logmoxregion[0,0,0] )) ; } 

{ 

Plot3 := CreatePlotClogplot') ; 
Plot(Plot3, logmoxregion); 

> 

{ Plot(Plot3, line, Polyvalue(a, line)); } 

sigma := Sqrt(Abs(-1.0 / (2.0 * a[2,0,0]))); 
{ SetTitle( Plotl, str (2*sigma) ) ; } 

WriteLn( 'gaussian.f it 2*sigma= ' , str (2*sigma) ) ; 

SetPlotStyle ( ps_Line ) ; 
SetPlotColor( col_Blue ) ; 

PlotXY(Plotl, linelOO, Exp(Polyvalue(a, linelOO))); 
SetPlotColor( col_Red ) ; 

PlotXY(Plotl, line, Exp(Polyvalue(a, line))); 
{SetPlotColor( col.Blue ) ; 

PlotXY(Plot2, linelOO, Polyvalue(a, linelOO));} 
end; 

{ Take moments } 

sumO := 0; suml := 0; sum2 := 0; 

for i :=0 to iylimit do 

begin 

sumO : = sumO + mox [0 , i , 0] ; 
suml := suml + mox[0,i,0] * i; 
sum2 : = sum2 + mox [0 , i , 0] * i * i ; 
end; 

WriteLn('y: avg=', suml/sumO, ' sigma^, 
Sqrt (sum2/sum0 - suml*suml/ (sum0*sum0) ) , ' sum=', sumO) ; 

{ WriteLn(textwin, 'max at y= ; , whereymax) ; } 



SetTitle( Plotl, StrCat (StrCat (StrCat (StrCat ( 'y : Colplot , max= ) , Str(micp)), 
StrCat( ; at ', Str (whereymax) ) ), ' sigma*2= '), Str (2*sigma) ) ); 

signal := micp; 
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cts_per_ms := total_counts_x / exposure_time ; 

WriteLn( 'blocker^ ,x, ' raicp=' , signal, ' total_cts_x= ; ' , total_counts_x, 
' cts_per_ms = ', cts_per_ms) ; 



data_x[index+l,0,0] := x; 
data_y[index+l,0,0] := signal; 
data_x [index, 0,0] :=x; 
data_y [index, 0,0] := signal; 



if do_dataplot then 
begin 

ClearPlot(Plot2) ; 

data_x_temp := data_x[0. . (index) ,0,0] ; 
data_y_temp := data_y[0. . (index) ,0,0] ; 

{ WriteLn(data_x[0. . index+1,0,0] , data.y [0 . . index+1 ,0,0] ) ; } 

PlotXY(Plot2, data_x_temp, data_y_temp) ; 

end; 

WriteLn( , frame= : ' , index, ' x=',x, ' micp=' , micp, ' mirp=', mirp, ' tot_x=' , total_co 
' tot_y', total_counts_y) ; 

{ ***** END OF FLUORESCENCE ***** } 

{ ***** TRIGGER RAMP AND COLLECT DATA ***** } 

if do_microlaser then 
begin 

SelectPort(2) ; 

Transmit ('A 80; { OUTPUT 1 -> HIGH } 
Transmit ('A 0'); { OUTPUT 1 -> LOW } 



SelectPort(l) ; 
RxFlush; 

Transmit ( 'CR' ) ; { SR400 Counter Reset } 

Transmit (' FB ' ) ; { Start with continuous data output } 

for i := 1 to n_periods do 
begin 
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dataline := RxString; 

data_array [i , index, 0] := val (dataline) ; 
{ WriteLn(data_array [i , index, 0] ) ; } 
end; 

data_max := 0; data_bkgnd := 0; 
data_bkgnd_sample : = 10; 

for i :=1 to n_periods do 
begin 

{ WriteLnCi, index=', i, ' ', index); } 
{ WriteLn(i, ' ', data_array) ; } 
if data_array [i , index, 0] > data_max then 
begin 

data_max := data_array [i, index, 0] ; 
where_max : = i ; 
end; 
end; 

data_bkgnd : = MeanOf (data_array [1 . . data_bkgnd_ sample , index , 0] ) ; 

WriteLn( 'data_max=' , data_max, ' at ' , where_max, '. data_bkgnd=' , 
data_bkgnd, ' signal=' , data_max-data_bkgnd) ; 

if (do_blackbody) then 

WriteLn(textwin, whereymax, ' ' , data_max - data_bkgnd) 
else 

WriteLn(textwin, index, ' ' , x, ' ', cts_per_ms, ' ', data_max - data_bkgnd) ; 
{ WriteLn(textwin, index, ' ', micp, } 

Clear(plot2) ; 

Plot (plot2 , data_array [1 . . n_periods , index , 0] ) ; 

if Save (data_array , dataf ilename) <> file_0k 
then Halt ('Save error'); 

x := x + x_increment; 
index : = index + 1 ; 

end; { end of 'if do_microlaser ' } 

until x > x_stop; 

{ *************** END OF MAIN LOOP ************** } 
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CloseSerial(l) ; CloseSerial (2) ; 



end 



C.2 Cavity locking experiment 



{ microlaser_experiment . v } 
{ C. Fang-Yen } 

{ Variables } 
var 

Zl ; plotimage ; bgimage ; img; 

PlotO ; Plotl ; Plot2 ; Plot3 ; expo sure .time ; 

imagename; backgroundname ; scatteringname ; scatteringimage; 

rowplotname ; colplotname ; 

do_dataplot; dataplotname ; crop_colplot ; 

i; index; x; { scanned voltage, etc. } 

x_start; x_stop; x_increment; x_use; 

micp; mirp; mox; moy; 

writedata; textwin; textwinname; 

do_serial; total_counts_x; total_counts_y ; num_pixels; 

wherexmax; whereymax; iylimit; ixlimit; sumO; suml; sum2; 

rezero; rezero_sample ; moymox; y_third_max_minus ; y_third_max_plus ; 

ii; do_gaussian_f it ; closest; fitmin; fitmax; 

sigma; logmox; logmoxregion; a; line; linelOO; 

test; absmox; data_x; data_y; data_x_temp ; data_y_temp ; 

doing_background_now; signal; previous_signal ; ratio; serial_error ; 

serial_error_msg; dataline; current_position; background_plot ; 

n_periods; data_array; data_array2; n_frames; data_max; data_bkgnd; 

where_max; 

data_bkgnd_sample ; v_t_553_on; v_t_791_on; v_t_553_off; v_t_791_off; 
use_background; plotimage_no553 ; dataf ilename ; do_blackbody ; do_microlaser ; 
save_images ; image_array; image_array_f ilename ; log_f ilename ; do_blocker; 
lock_time; do_553_ttl; cts_per_ms; 

procedure ttl_outputs (outl , out2, out3) ; 
begin 

SelectPort(2) ; 

Transmit ('A '+ Str (8* (1-outl) + 16*(l-out2) + 32* (l-out3) ) ) ; 
end; 
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procedure lock_cavity; 
begin 

ttl_outputs(0 , 0, 0); { turn off AOM after cavity } 
ttl_outputs (1 , 1, 0); { turn on 791 probe; no cav lock disable} 
end; 

procedure unlock_cavity ; 
begin 

ttl_outputs (0 , 0, 0); { turn off 791 probe; cav lock disable } 

ttl_outputs (0 , 0, 1); { turn on AOM after cavity } 

end; 

procedure take_picture(roi_index; var img) ; 
begin 

{ Set Gain Index 3 } 
pvcSetGain( 3 ) ; 
exposure_time := 300; 
pvcSetExpTime ( exposure_time) ; 
roi_index := 0; { override } 

if(roi_index = 0) then 

img := pvcSequence( 1,386,228,423,257,1,1 ) ; {small} 
if(roi_index = 1) then 

img := pvcSequence( 1,304,169,501,323,1,1 ) ; {medium} 
if(roi_index = 2) then 

img := pvcSequence( 1,0,0,767,511,1,1 ); { full CCD } 

if(roi_index = 3) then 

img := pvcSequence( 1,352,236,454,270,1,1 ); { small for 
slit} 

end; 
begin 

-[ ******************* INITIALIZATION **************** } 

datafilename := ' data.tif; { ' data.tif = generic } 
image_array_f ilename := ' image .tif ' ; 
log_filename := ' log.txt'; 

do_blackbody := False; { Record location of blackbody peak in col plot } 
use_background := False; { False = use scattering subtraction } 
do_microlaser := False; { False = skip cavity scan } 
save_images := False; { save fluorescence images } 
do_blocker := True; 

imagename : = ' image ' ; 
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rowplotname := 'rowplot'; 

colplotname := 'colplot'; 

dataplotname := 'dataplot ' ; 

crop_colplot := 0; 

backgroundname : = ' background ' ; 

scatteringname := 'scattering'; 

index := 0; 

x_use := True; 

x_start := 20; 

x_stop := 60; 

x_increment := 1; 

x := x_start; 

writedata := True; 

textwinname : = ' output ' ; 

do_serial := True; 

rezero := False; 

rezero_sample := 10; 

do_gaussian_f it := False; 

data.x := MakeLinear (0 . , . 0, 1000) ; 

data.y := MakeLinear (0 . , . 0, 1000) ; 

do_dataplot := False; 

doing_background_now := 0; 

v_t_553_on := 0.605; { V_t 553nm on voltage } 
v_t_791_on := 1.0; { V_t 791nm on voltage } 
v_t_553_off := -0.1; { V_t 553nm off voltage } 
v_t_791_off := -0.1; { V_t 791nm off voltage } 

lock.time := 1000; 



if (writedata) then 
begin 

if EditorExists (textwinname) then 
GetEditor (textwinname , textwin) 
else textwin := CreateEditor ( 'output' ); 

WriteLn(textwin, ' ' ) ; 

end; 

{ ******* INITIALIZE SERIAL PORTS ****** } 

SetTxEnd( chr(13) ) ; 
SetRxEnd( chr(13) ) ; 
SetRxTimeout(lOOO) ; 

WriteLnC ') ; 
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OpenSeriaK 1 , 19200,8,NoParity, 1 ) ; { SR400 photon counter A } 
serial_error := SerialError; 

if serial_error <> then 
Halt( 'Failed to open serial connection 1: error ', serial_error) 
else WriteLn( 'Serial connection 1 ok'); 

SelectPort(l) ; 
TxFlush ; 

Transmit ( chr(13)) ; 

Transmit ( chr(13)) ; 
TxFlush ; 
RxFlush; 

Transmit (' pi 1'); 

dataline := RxString; 
WriteLn('SR400A: P0RT1=' , dataline); 

Transmit ( 'np' ) ; 

dataline := RxString; 
n_periods := val (dataline) ; 
WriteLn('SR400A: NPERI0DS=', n.periods) ; 

{ 

CloseSerial(l) ; 

Halt ( " ) ; 

} 

if x_increment <> then 

n_frames := Floor ( (x_stop - x_start)/x_increment)+l 
else 

n_frames := 100; 

data_array := CreateArray (single, n_frames, 1); 
data_array2 := CreateArray (single, n_frames, 1); 

OpenSeriaK 2,9600,8,NoParity, 1 ) ; { IMS stepper driver } 
serial_error := SerialError; 

if serial_error <> then 
Halt( 'Failed to open serial connection 2: error ', serial_error ) 
else WriteLn( 'Serial connection 2 ok'); 
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{ ****** INITIALIZE BLOCKER ****** } 

SelectPort(2) ; 

Transmit ( ' ' ) ; 

Delay(lOO) ; 
RxFlush; 

{ WriteLn( 'RxWaiting^ ,RxWaiting) ; } 
Transmit ( ' Z ' ) ; 
Delay (100) ; 
dataline := RxString; 

WriteLn ( 'Current position dataline); 
dataline := StrParse (dataline , 'Z'); 
current_position := Val (StrStrip (dataline) ) ; 
WriteLn( 'Current position Z=', current_position) ; 
Transmit (' X' ) ; 
dataline := RxString; 
WriteLn (dataline) ; 
dataline := RxString; 
WriteLn (dataline) ; 

{ ************* TAKE 1 FRAME, SET UP IMAGE AND PLOT WINDOWS ************ } 
unlock_cavity ; { unlock cavity } 
take_picture(0, plotimage) ; 
lock_cavity; { re-lock cavity} 
Show( plotimage) ; 

num_pixels := GetXSize (plotimage) * GetYSize (plotimage) ; 
WriteLn ( 'n_frames=' , n_f rames) ; 
if save .images then 

image_array := CreateArray (GetXSize (plotimage) , 
GetYSize (plotimage) , n_f rames) ; 

WriteLn (FileExpand( 'background.tif ') ) ; 
{ Halt; } 

if use_background then 
if not ImageExists(backgroundname) then 
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begin 

if Query ('No background: use this? ' )=id_Yes then 
begin 

Save (plot image, 'background.tif') ; 
Open ( 'background.tif ' , background_plot) ; 
Show (background_plot , ' background ' ) ; 
Delete (plotimage) ; 
end; 
Halt; 
end; 

{ 

if use_background = False then 

if not ImageExists(scatteringname) then Halt('No scattering file'); 

} 

{ Get the active frame number } 

Zl := GetActiveFrame ( plotimage ) ; 

{ Create a plot window } 
if PlotExists (rowplotname) then 
begin 

GetPlot (rowplotname , PlotO) ; 
{ Halt ('ok'); } 
end 

else 
begin 

PlotO := CreatePlot( rowplotname ); 
SetTitle( PlotO, 'Row Plot' ) ; 
SetSubTitle( PlotO, GetName( plotimage ) ) ; 
SetXLabeK PlotO, 'X Position' ) ; 
SetYLabeK PlotO, 'Mean Intensity' ) ; 
SetGrid( PlotO, gs_Both ) ; 
end; 

{ Row plot } 

SetPlotColor( col_Black ) ; 

Plot( PlotO, MeanOverY( plotimage [..,.., Zl] ) ) ; 

{ Create data plot window } 
if PlotExists (dataplotname) then 
begin 

GetPlot (dataplotname , Plot2) ; 
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end 
else 

Plot2 := CreatePlot (dataplotname) ; 

{ Create a plot window } 
if PlotExists (colplotname) then 
begin 

GetPlot (colplotname, Plotl) ; 

{ Halt (' ok'); } 
end 
else begin 

Plotl := CreatePlot ( colplotname ) ; 

SetTitle( Plotl , 'Column Plot' ) ; 

SetSubTitle( Plotl,GetName( plotimage ) ) ; 

SetXLabeK Plotl, 'Y Position' ) ; 

SetYLabeK Plotl, 'Mean Intensity' ) ; 

SetGrid( Plotl ,gs_Both ) ; 
end; 

{ Col plot } 

SetPlotColor( col_Black ) ; 

Plot( Plotl ,MeanOverX( plotimage [..,.., Zl] ) ) ; 

{ *** INITIALIZE SR400 PARAMETERS *** } 
{ to be written } 

{ ******************* MAIN LOOP BEGINS ************************ } 
repeat 

if (do_blocker) then 
begin 

SelectPort(2) ; 

Transmit (StrCatCR ', Str(x))); { SET BLOCKER POSITION } 

Delay (1000. 0*Abs (current.position - x)/400.0); { Wait for it to move } 

current_position := x; 

end; 

{ Turn on 553nm probe, turn off 791nm pump} 
SelectPort(l) ; 

Transmit (StrCatC pi 1, ', Str (v_t_553_on) ) ) ; 
Transmit (StrCat ('pi 2, ', Str (v_t_791_of f ) ) ) ; 
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unlock_cavity ; {unlock} 



{ Acquire data } 
take_picture(0, plotimage) ; 
Show( plotimage ) ; 

lock_cavity; {relock + Turn off 553nm probe} 

{ Turn off 553nm probe } 
SelectPort(l) ; 

Transmit (StrCatO pi 1, Str (v_t_553_of f ) ) ) ; 

{ Transmit (StrCatO pi 2, ', Str (v_t_791_on) ) ) ; } 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 

{ plotimage [.., .. crop_colplot ,.. ] := mox[0,crop_colplot,0] ; } 

Get Image ( backgroundname , bgimage ) ; 
{ 

Getlmage( scatteringname , scatteringimage ); 
} 

{ Perform an image subtraction } 

if use_background then plotimage := plotimage - bgimage ; 

if not use_background then 
begin 

Delay (lock_time) ; 

unlock_cavity ; { unlock } 

take_picture(0, plotimage_no553) ; 

lock_cavity; { relock } 
{ Halt (' should be locked');} 

Delay (lock_time) ; 
{ plotimage := plotimage - plotimage_no553 - scatteringimage;} 

plotimage := plotimage - plotimage_no553 ; 

end; 

{ Rezero: subtract avg of first 10 pixels in col plot } 
if rezero then begin 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 
moymox : = MeanOverY (mox [ . . , . . 10 , . . ] ) ; 
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plotimage := plotimage - moymox [0 , , 0] ; 

{ mox[0,crop_colplot,0] ; } 

end; 

if (save_images) then 
begin 

image_array [..,.., index] := plotimage; { store plotimage in array } 
if Save (image_array , image_array_f ilename) <> file_0k 
then Halt ('Save error'); 

end; 

{ Row plot } 

ClearPlot(PlotO) ; 

SetPlotColor( col_Black ) ; 

moy := MeanOverY( plotimage [..,.. ,Z1] ); 

Plot( Plot0,moy ) ; 
mirp := MaxOf (moy) ; 

wherexmax := 0; 

ixlimit := GetXSize (moy) -1 ; 

for i := to ixlimit do 

begin 

if moy [i, 0,0] = mirp 

then wherexmax := i; 
{ WriteLn(str(i)) ; } 
end; 

{ 

sumO := 0; suml := 0; sum2 := 0; 

for i :=0 to ixlimit do 

begin 

sumO := sumO + moy [i, 0,0] ; 
suml := suml + moy [i, 0,0] * i; 
sum2 := sum2 + moy [i, 0,0] * i * i; 
end; 

WriteLn('x: avg=', suml/sumO, ' sigma=', 
Sqrt (sum2/sum0 - suml*suml/ (sum0*sum0) ) , ' sum=', sumO) ; 

WriteLn(x, ' ' , sumO) ; 
} 

{ WriteLn(textwin, 'max at x=', wherexmax); } 

{ WriteLn(textwin, MaxOf (moy [0 . . 10] ) ) ; } 
total_counts_y := MeanOf (moy) * num_pixels; 
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SetTitleC PlotO, StrCat (StrCat ( 'Row Plot, max=', Str(mirp)), 
StrCat (' at Str (wherexmax) ) ) ) ; 

{ Col plot } 

ClearPlot(Plotl) ; 

SetPlotColor( col_Black ) ; 

mox := MeanOverX( plotimage [..,.. ,Z1] ); 

{ WriteLn(textwin, GetYSize(mox)) ; } 

mox [ . . ,0 . . crop_colplot , . . ] : = mox [ . . , crop_colplot , . . ] ; 
{ if crop_colplot>0 then 

for i := 1 to crop_colplot-l do 
begin 
WriteLn(i) ; 

mox[i] := mox [crop_colplot] ; 
end; } 

Plot( Plotl,mox ) ; 
micp := MaxOf (mox) ; 

total_counts_x := MeanOf (mox) * num_pixels; 

whereymax := 0; 

iylimit := GetYSize (mox) -1 ; 

for i := to iylimit do 

begin 

if mox[0,i,0] = micp 

then whereymax := i; 
{ WriteLn(str(i)) ; } 
end; 

{ Find 1/3 max points } 

closest := 10000; 
for ii := whereymax to iylimit do 
begin 

if Abs(mox[0,ii,0] - micp/3.0) < closest then 
begin 

closest := Abs(mox[0,ii,0] - micp/3.0); 
y_third_max_plus := ii; 
end; 
end; 

closest := 10000; 
for ii:=0 to whereymax do 
begin 

if Abs(mox[0,ii,0] - micp/3.0) < closest then 
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begin 

closest := Abs(mox[0,ii,0] - micp/3.0); 

y_third_max_minus := ii; 

end; 
end; 
{ 

WriteLn( 'yjialf _max_plus str (y_third_max_plus) ) ; 
WriteLn( ; y_half _max_minus str (y_third_max_minus) ) ; 
} 

{ Gaussian fitting } 
sigma :=0; 

if do_gaussian_f it then begin 

absmox := sqrt(sqr(mox))+0.01; 
logmox := Ln (absmox) ; 
fitmin := y_third_max_minus ; 
fitmax := y_third_max_plus ; 

line := MakeLinear (fitmin -0.0 , fitmax - 0.0, fitmax - fitmin + 1) ; 
linelOO := MakeLinear (0, iylimit+0 . 0, (iylimit+1 . 0) ) ; 
a := PolyFit (line , logmox [0, fitmin. .fitmax, 0] ,2) ; 

{ WriteLn(a) ; } 
{ 

WriteLn(str (line) ) ; 

WriteLn(str(logmox[0, fitmin. .f itmax,0] )) ; 
WriteLn(str(a[0,0,0] )) ; 
WriteLn(str(a[l,0,0] )) ; 
WriteLn(str(a[2,0,0] )) ; 
WriteLnO '); } 

{ logmoxregion := logmox [0, fitmin. .fitmax, 0] ; 
WriteLn(logmoxregion) ; 
WriteLn(str(logmoxregion[0,0,0] )) ; } 

{ 

Plot3 := CreatePlotClogplot') ; 
Plot(Plot3, logmoxregion); 

} 

{ Plot(Plot3, line, Polyvalue(a, line)); } 

sigma := Sqrt(Abs(-l .0 / (2.0 * a[2,0,0]))); 
{ SetTitle( Plotl, str (2*sigma) ) ; } 
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WriteLn( 'gaussian_f it 2*sigma= ' , str (2*sigma) ) ; 

SetPlotStyle ( ps_Line ) ; 
SetPlotColor( col_Blue ) ; 

PlotXY(Plotl, linelOO, Exp(Polyvalue(a, linelOO))); 
SetPlotColor( col_Red ) ; 

PlotXY(Plotl, line, Exp(Polyvalue(a, line))); 
{SetPlotColor( col.Blue ) ; 

PlotXY(Plot2, linelOO, Polyvalue(a, linelOO));} 
end; 

{ Take moments } 
{ 

sumO := 0; suml := 0; sum2 := 0; 

for i :=0 to iylimit do 

begin 

sumO : = sumO + mox [0 , i , 0] ; 
suml := suml + mox[0,i,0] * i; 
sum2 : = sum2 + mox [0 , i , 0] * i * i ; 
end; 

WriteLn('y: avg=', suml/sumO, ' sigma^, 
Sqrt (sum2/sum0 - suml*suml/ (sum0*sum0) ) , ' sum=', sumO) ; 

> 

{ WriteLn(textwin, 'max at y=', whereymax) ; } 

SetTitle( Plotl, StrCat (StrCat (StrCat (StrCat ( 'y : Colplot , max=', Str(micp)), 
StrCat (' at Str (whereymax) ) ), ' sigma*2= '), Str (2*sigma) ) ); 

signal := micp; 

cts_per_ms := total_counts_x / exposure_time ; 

WriteLn( 'blocker^ ,x, ' raicp=' , signal, ' total_cts_x= ; ' , total_counts_x, 
'cts.per.ms = cts_per_ms) ; 



if do_dataplot then 
begin 

data_x[index+l,0,0] := x; 
data_y[index+l,0,0] := signal; 
data_x [index, 0,0] :=x; 
data_y [index, 0,0] := signal; 

ClearPlot(Plot2) ; 
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data_x_temp := data_x[0. . (index) ,0,0] ; 
data_y_temp := data_y[0. . (index) ,0,0] ; 



{ WriteLn(data_x[0. . index+1,0,0] , data.y [0 . . index+1 ,0,0] ) ; } 

PlotXY(Plot2, data_x_temp, data_y_temp) ; 

end; 

WriteLn('frame=' , index, ' x=',x, ' micp=' , micp, ' mirp=' , mirp, ' tot_x=' 
' tot_y', total_counts_y) ; 

{ ***** END OF FLUORESCENCE ***** } 

{ ****** MICROLASER I ******* } 

unlock_cavity ; { unlock } 

{ Delay (100); } { Wait a bit } 

SelectPort(l) ; 

Transmit (StrCatO pi 2, ', Str (v_t_791_on) ) ) ; { turn on pump AFTER unlock } 

SelectPort(l) ; 
RxFlush; 

{ Transmit OCR') ; }{ SR400 Counter Reset } 
{ Transmit ( 'FB' ) ; }{ Start with continuous data output } 
Delay(lOO) ; 

Transmit (' QB' ) ; { Read most recent complete data point } 

dataline := RxString; 

data_array [index, 0] := val (dataline) ; 

data_bkgnd := 0; { Enter correct value here } 

if (do_blackbody) then 

WriteLn(textwin, whereymax, ' ', data_array [index, 0] - data_bkgnd) 
else 

WriteLn(textwin, index, ' ' , x, ' ' , micp, ' ' , 
total_counts_x/exposure_time, ' ', 
Str (data_array [index, 0] - data_bkgnd) ) ; 

{ Clear (plot2) ; 

Plot(plot2, data_array [1 .. index, 0] ) ; } 

{ if Save(data_array, dataf ilename) <> file_0k 
then Halt ('Save error'); } 
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lock_cavity; { relock } 
Delay (lock_time) ; 

unlock_cavity ; { unlock, pump still on! } 
-[ *********** MICROLASER II ******** } 

SelectPort(l) ; 
RxFlush; 

{ Transmit OCR') ; }{ SR400 Counter Reset } 
{ Transmit ( 'FB' ) ; }{ Start with continuous data output } 
Delay(lOO) ; 

Transmit (' QB' ) ; { Read most recent complete data point } 

dataline := RxString; 

data_array2 [index, 0] := val (dataline) ; 

data_bkgnd := 0; { Enter correct value here } 

if (do_blackbody) then 

WriteLn(textwin, whereymax, ' ', data_array [index, 0] - data_bkgnd) 
else 

WriteLn(textwin, index, ' ' , x, ' ' , micp, ' ' , 
total_counts_x/exposure_time, ' ', 
dat a_ array 2 [index, 0] - data_bkgnd) ; 

lock_cavity ; 
Delay (lock_time) ; 

{ 

Clear(plot2) ; 

Plot(plot2, data_array [1 .. index, 0] ) ; } 

{ if Save(data_array2, dataf ilename) <> file_0k 
then Halt ('Save error'); 

} 

x := x + x_increment; 
index : = index + 1 ; 

until x > x_stop; 

{ *************** END OF MAIN LOOP ************** } 
CloseSerial(l) ; CloseSerial(2) ; 
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end 
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Appendix D 

Standing-wave versus 
traveling-wave coupling 



We now show that the peak atom-cavity coupling g in the traveling-wave case is equal 
to one-half the peak coupling of the standing- wave case. 

Consider one atom in a cavity, located along the cavity's axis at z — z(t) where 



z(t) = Zq + V t 



(D.l) 



We can write the operator for the total electric field seen by the atom: 
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Here V = -kLwq/A is the mode volume, with L the cavity mirror spacing. 

Writing out z(t), we can see the splitting of the two cavity frequencies ui ± kv: 
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Suppose v — 0. Then 
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This is the familiar standing wave field operator. 

For tJ^Owe can set z = by redefining the time origin. (This works as long as 
kvt int > 2n) 
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Suppose that only one these terms is resonant. Since kg = ft ■ E, we have shown 
explicitly that the traveling-wave case has a maximum coupling one-half that of the 
standing-wave case. 
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